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1.  INTRODUCTION 


1.1.  Object  and  Scope 

The  main  objective  of  this  study  is  to  develop  a  set  of  numerical 
methods  suitable  for  investigating  the  load-deflection  and  bifurcation 
characteristics  of  structures  for  which  significant  nonlinear  behavior  is 
possible.  The  methods  are  applicable  to  a  wide  variety  of  structures,  but 
will  be  examined  in  detail  only  with  reference  to  one  .  f  the  simplest  types 
of  structures  possessing  the  necessary  complications  i.n  behavior  -  the 
planar  arch  under  a  concentrated  load. 

The  term  "planar",  as  used  in  this  study,  refers  to  the  configu¬ 
ration  of  the  arch  during  the  initial  stages  of  loading  (often  called  the 
prebuckling  configuration) .  Both  in-plane  and  out-of-plane  buckling  be¬ 
havior  of  the  planar  configuration  are  examined.  Although  it  would  be 
possible  to  include  the  effect  of  certain  nonlinear  stress-strain  laws, 
the  nonlinear  behavior  examined  in  this  study  is  geometrical  in  nature 
and  results  from  large  displacements  (arising  from  large  rotations  but 
small  strains). 

The  numerical  methods  developed  here  are  capable  of  determining 
limit  points  on  the  load-deflection  curve  (see  Fig.  2,  points  A  and  B) , 
as  well  as  finding  bifurcation  points  and  subsequently  tracing  the  buck¬ 
led  configuration.  The  numerical  results  given  in  Chapter  5  illustrate 
these  capabilities  in  problems  of  considerable  technical  interest. 

1.2.  General  Remarks  and  Observations 

From  the  errliest  work  on  the  buckling  of  cylindrical  shells, 
it  has  been  nottd  that  experimentally  determined  buckling  loads  of  various 
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In  view  of  this  wide  variety  of  possible  behavior  of  structural 
members  a  consideration  of  postbuckling  behavior  is  an  essential  part  of 
the  analysis  of  a  given  structure  which  exhibits  a  buckling  phenomenon. 

1.3.  Background 

As  mentioned  above,  the  numerical  methods  developed  in  the  pre¬ 
sent  study  are  applied  to  the  simplest  structures  which  exhibit  the  non¬ 
linear  behavior  necessary  to  provide  an  adequate  test  of  the  methods.  The 
mathematical  model  of  the  structures  studied  here  is  given  by  Love  (1927) 
for  the  equilibrium  forms  of  thin  rods.  According  to  Love,  Clebsch  (1862) 
and  Kirchloff  (1859)  arrived  independently  at  the  equilibrium  equations. 
The  geometrical  relationships  are  attributed  to  Routh  (1905) ,  and  Clebsch 
(1862)  is  given  credit  for  the  moment -curvature  relationships.  These 
equations  presented  by  Love  are  applicable  to  the  three-dimensional  be¬ 
havior  of  thin,  linearly  elastic  rods  with  inextensional  centerlines,  al¬ 
though  an  indication  is  given  by  Love  of  the  necessary  modification  for  an 
extensional  centerline.  Vlasov  (1959)  indicates,  that  as  a  first  approxi¬ 
mation,  the  effect  of  warping  restraint  on  the  behavior  of  curved  beams 
may  be  introduced  by  using  the  corresponding  relationship  between  torque 
and  rate  of  twist  for  a  straight  rod.  In  Chapter  5,  results  are  presented 

P/\  v»  ►Ka.  ******  UnoUl  a  (•  ovnlioe  ttUavq  f  Vn  aPPaa**  a  4-  av^aao  J  a«  a  f  ►V.a 

i.V  b  WltW  J.14  £>  V 1  Vtt.  W4IWU  n»»WbW  Wt»W  Vbb’^Vt.  Vb  V  b 

centerline  is  included  and  for  the  lateral  buckling  of  an  I-beam  where 
warping  restraint  is  considered. 

The  oldest  analysis  of  buckling,  Eller's  work  on  a  perfect 
elastic  column,  (see  Timoshenko  and  Gere  (1961))  included  a  postbuckling 
analysis.  ‘'However.,  the  perfect  column  is  one  case  in  which  the  behavior 
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results  for  various  rise-to-span  ratios.  The  mathematical  model  assumed 
an  inextensional  centerline.  It  is  not  clear  whether  or  not  extension 
of  the  centerline  would  complicate  this  computational  method,  which  in¬ 
volved  elliptic  integrals. 

1.4.  Outline  of  the  Method  of  Analyses 

In  this  study  a  set  of  numerical  techniques  is  developed  for 
improving  an  approximation  to  a  bifurcation  point  on  the  load-deflection 
curve.  One  method  permits  a  direct  computation  of  an  approximate  eigen¬ 
vector  which  is  then  improved  simultaneously  with  the  prebuckling  config¬ 
uration. 

The  technique  requires  a  solution  of  a  set  of  nonlinear  equations 
which  indicate  how  the  prebuckling  configuration  (including  the  loading) 
must  be  modified  in  order  to  reach  the  bifurcation  point.  This  part  of 
the  solution  is  treated  in  Chapter  2  in  a  mathematical  fashion  and  in 
Chapter  4  for  a  specific  physical  problem.  The  nonlinear  equations  are 
developed  with  reference  to  the  general  eigenvalue  problem  A  X  =  X  B  X 
and  are  solved  by  a  modification  of  the  Newton-  jhson  method. 

As  indicated,  the  solution  process  predicts  how  the  prebuckling 
configuration  must  be  changed  to  reach  a  bifurcation  point.  The  process 
of  modifying  the  prebuckling  configuration  is  examined  in  Chapter  3.  The 
standard  NewtGn-Raphson  procedure  may  be  used  except  when  the  prebuckling 
configuration  is  near  a  bifurcation  point.  As  noted  by  Thurston  (1969), 
the  equations  specifying  the  linear  changes  in  the  prebuckling  configura¬ 
tion  become  singular  at  bifurcation  points.  A  method  proposed  in  this 
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study  actually  makes  use  of  this  fact  to  arrive  at  an  improved  -prebuckling 
configuration  and  a  better  estimate  of  the  eigenvector  in  a  rapidly  con¬ 
vergent  computation. 

1.5.  Nomenclature 


The  symbols  used  in  this  study  are  defined  in  the  text  when  they 
first  appear.  For  convenient  reference,  tne  more  important  symbols  are 
summarized  here  in  alphabetical  order.  Some  symbols  are  assigned  more  than 
one  meaning:  however,  in  the  context  of  their  use  there  are  no  ambiguities. 


a 

A,  B,  C 

b 


C,  C,  1),  D 
det(x) 


El. 

i 

Ep  -5EI 


EC 


W 


H 


I 


1 


I 


2 


I 


3 


radius  of  undef  rmed  circular  arch 

general  liaearlized  operators,  may  be  matrices 
differential  or  Integral  operators 

constant  vector 

coefficient  matrices  of  linear  algebraic  equations 
determinant  of  x 

deflection  components  at  concentrated  load,  in 
global  coordinates  i  =  1,  2,  3 

scalar  error  term 

flexural  rigidities  (includes  St.-Venant  tor¬ 
sional  rigidity),  i  «  1,  2,  3 

the  Itn  configuration  and  its  corresponding 
increment  in  the  Newton-Raphson  procedure 

warping  rigidity 


rise  of  undeformed  arch 

for  a  planar  member,  moment  inertia  about  an 
axis  perpendicular  to  the  plane 

for  a  planar  member,  moment  of  inertia  about 
an  axis  in  the  plane 


corresponds  to  J,  the  St.-Venant  torsion  constant 


2.  PROCEDURE  FOR  FINDING  BIFURCATIONS 


2.1.  General 

A  study  of  pcstbuckling  behavior  requires  at  least  two  items  of 
information.  These  are  the  buckling  load,  along  with  the  corresponding 
configuration  just  prior  to  buckling,  and  the  eigenvector,  which  gives  an 
initial  estimate  of  the  postbuckling  path.  In  the  following  sections 
theoretical  considerations  are  presented  which  lead  to  the  development  of 
a  set  of  efficient  numerical  methods  for  treating  bifurcations  from  a 
nonlinear  prebuckling  state.  Detailed  descriptions  of  the  numerical  pro¬ 
cedures  are  reserved  for  Chapters  3  and  4. 

2.2.  Bifurcation  as  an  Eigenvalue  Problem 


The  eigenvalue  problems  to  be  treated  here  are  assumed  to  be 
described  by 


A  X  =  X  B  X 


(2.1) 


and  appropriate  boundary  conditions  where  necessary.  The  quantities  A  and 
B  may  be  matrices,  differential,  or  integral  operators;  X  is  the  eigen¬ 
value  and  X  the  eigenvector.  The  operators  A  and  B  refer  to  the  prebuck¬ 
ling  configuration  and  are  in  general  dependent  on  the  eigenvalue  X  but 
not  on  the  eigenvector  X.  It  is  assumed  that  the  dependence  of  A  and  B  on 
X  is  known,  at  least  implicitly. 

The  discrete  (algebraic)  eigenvalue  problem  may  be  represented 
by  Eq.  (2.1)  when  A  and  B  are  interpreted  as  matrices.  One  technique 
that  has  been  used  to  solve  this  type  of  problem  is  to  increment  the  trial 


$ 


a 


3 


5 

t.  M 

i3k 


A,  A,  A 
*  cr 

* 


noi.-dimensionalized  buckling  load  (out-of-plane), 


non-dimensionalized  buckling  load  (in-plane) 

8  =  Pa2/EI1 

increment  operator 
alternating  tensor 

strain  of  centerline 

eigenvalues 

used  to  denote  eigenvector  quantities 
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eigenvalue  A  (which  in  general  implies  changing  a  and  B)  and  at  each 
value  of  A  to  compute  the  determinant  of  (A  -  AB) .  This  procedure  was 
used  by  Leicester  (1968)  and  in  essence  is  an  extension  of  the  so-called 
Holzer  method,  Holzer  (1921).  A  change  of  sign  of  this  determinant 
between  successive  values  of  the  trial  eigenvalue  indicates  an  eigenvalue 
falling  in  that  range.  Interpolation  may  be  used  to  find  the  value  of  A 
for  which  det  (A  -  AB)  is  zero.  At  this  stage,  the  eigenvector  may  be 
generated  in  the  conventional  manner  by  setting  one  of  the  components  of 
X  to  unity  (say  X^)  and  solving  for  the  other  components  on  this  basis. 

It  may  be  appropriate  to  mention  that  det  (A  -  AB)  equal  to  zero  does 
not  necessarily  imply  bifurcation.  It  may  mean  that  there  is  a  limit 
point  on  A,  and  some  other  quantity  should  be  incremented. 

The  linearized  equation  governing  the  local  behavior  of  the 
branch  of  the  equilibrium  curve  corresponding  to  the  prebuckling  config¬ 
uration  is  of  the  form  (A  -  AB)  Y  =  b.  It  is  then  evident  from  Eq.  (2.1) 
that  an  impending  singularity  of  (A  -  AB)  will  cause  numerical  difficulties 
associated  with  changing  the  prebuckling  configuration  in  the  vicinity  of 
a  bifurcation  point.  That  is,  changes  in  A,  B,  and  Y  will  not  be  accurate. 
This  has  been  noted  previously  by  Thurston  (1969) ,  who  presented  a  compu¬ 
tational  device  for  the  solution  in  that  case.  This  same  phenomenon  has 
been  encountered  in  this  study  and  the  means  of  computation  which  has  been 
devised  is  introduced  in  the  next  section.  It  will  be  seen  to  be  less  in¬ 
volved  than  that  presented  by  Thurston. 

The  continuous  eigenvalue  problem  may  be  solved  in  a  manner 
similar  to  the  discrete  problem.  In  this  case,  however,  it  is  not  det 
(A  -  AB)  which  is  examined  but  rather  the  determinant  corresponding  to 
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satisfaction  of  the  boundary  conditions.  This  technique  has  been  used 
by  Cohen  (1965),  Kalnins  (1964)  and  Zarghamee  and  Robinson  (1967).  As 
with  the  discrete  problem,  there  may  be  numerical  difficulties  in  deter¬ 
mining  accurate  changes  in  the  prebuckling  configuration  near  bifurca¬ 
tion  points. 

2.3.  A  New  Solution  Technique 

An  essential  characteristic  of  the  technique  presented  here  is 
the  simultaneous  Improvement  of  the  bifurcation  point  (load  and  configura¬ 
tion)  and  the  eigenvector  by  a  n-ocess  involving  the  interaction  between 
the  two. 

If  the  A,  B,  and  A  corresponding  to  a  particular  prebuckling 
configuration  and  an  approximate  eigenvector  are  substituted  into  Eq. 

(2.1),  then 

(  A  X  -  A  B  X  (2.2) 

th 

where  the  superscript  j  indicates  the  j  ‘  approximation  and  R  is  a 
residual.  The  object  then  is  to  remove  the  residual  from  Eq.  (2.2).  In 
the  usual  eigenvalue  problem,  A  is  not  treated  as  an  unknown  of  the  same 
type  as  X.  However,  the  method  proposed  here  considers  A  B  X  as  a  non¬ 
linear  term.  This  suggests  that  some  nodification  of  the  well-known 
Newton-Raphson  procedure  may  be  applicable  here.  Use  of  the  standard 
Newton-Raphson  technique  has  been  discussed  by  Kalnins  and  Lestingi 
(1967),  Leicester  (1968)  and  West  and  Robinson  (1969).  In  order  to 
extend  the  Newton-Raphson  technique  to  bifurcation  problems,  it  is 
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necessary  to  linearize  Eq.  (2.1)  about  some  known  configuration  (say  the 
j^).  in  essence,  Eq.  (2.1)  is  expanded  about  the  j ^  configuration  and 
only  the  linear  terms  are  kept. 

The  linearization  of  Eq.  (2,1)  yields 

(A6X  -  XB6X}^  =  {-6AX  +  6ABX  +  X6BX  -  R}^>  (2.3) 

Since  A  and  B  are  in  general  dependent  on  X,  the  linear  parts  of  the 
increments  of  A  and  B  may  be  formally  expressed  as 

I  < 

SA  =  |£i6X  ,  6B=|r!5X  (2.4) 

X=X(j) 

Substitution  of  Eqs.  (2.4)  into  Eq.  (2.3)  results  in 

(A6X  -  XB5X}(;5)  =  {5X(-  X  +  BX  +  X  X)  -  R}(j)  (2.5) 

Examination  of  Eq.  (2.3)  reveals  there  are  two  types  of  incremental  quan¬ 
tities  to  be  considered;  those  corresponding  to  changes  in  the  eigenvector 
6X  and  those  corresponding  to  changes  in  the  prebuckling  configuration  6X, 
6B,  and  6A,  From  Eq.  (2.4),  5A  and  5B  are  related  to  6X  so  that  in  fact, 
the  unknowns  are  dX  and  dA,  as  indicated  in  Eq.  (2.5). 

Once  the  quantities  ~,  ~  and  an  approximate  eigenvector  are 
computed,  the  solution  of  Eq.  (2.5)  may  proceed  as  follows.  Since  <5X  is 
an  unknown,  there  is  one  more  unknown  than  there  are  equations  to  solve, 
a  situation  that  does  not  arise  in  the  usual  Newton-Rapnson  technique. 

The  presence  of  ar.  extra  unknown  is  to  be  expected,  since  the  amplitude 


I  X=A^ 
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of  the  eigenvector  is  indeterminate.  The  arbitrariness  in  the  eigenvector 
is  removed  by  specifying  a  scalar  side  condition 

XT  B  <SX  =  0  (2.6) 

This  side  condition  (  or  its  integral  equivalent  when  appropriate)  allows 
a  solution  for  5X  and  6X  by  eliminating  the  possibility  cf  large  changes 
in  the  eigenvector  if  the  eigenvalue  and  approximate  eigenvector  are 
nearly  correct. 

If  the  computed  61  is  not  satisfactorily  small,  the  prebuckling 
configuration  is  not  one  corresponding  to  an  eigenvector  and  must  be  modi¬ 
fied.  The  magnitude  of  6A  dictates  how  the  procedure  continues.  In  es¬ 
sence,  this  method  predicts  approximately  how  A  and  the  prebuckling  con¬ 
figuration  should  be  changed  to  approach  a  bifurcation  point. 

For  the  above  solution  process,  it  has  been  implicitly  assumed 
that  the  quantities  could  be  computed.  From  Eq.  (2 . 5  >  it  appears 

that  these  quantities  could  be  obtained  by  computing  6A  and  6B  for  a  unit 
value  of  6A  (6A  =  1).  This  is  a  straightforward  application  of  the 
Newton-Raphson  procedure.  However,  as  mentioned  in  Chapter  1,  the  equa¬ 
tions  become  singular  at  bifurcation  points.  This  means  that  at  or  near 
bifurcation  points,  a  special  computational  davice  must  be  incorporated 
into  the  Hewton-Raphson  technique  in  order  to  compute  changes  in  the  pre¬ 
buckling  configuration  accurately.  This  special  computational  device  is 


discussed  in  the  next  section. 
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2.4.  Numerical  Treatment  of  the  Singular  Equations 

As  mentioned  above,  the  direct  procedure  for  changing  the  pre¬ 
buckling  configuration  is  bound  to  fail  at  or  near  the  bifurcation  point. 

The  difficulty  is  caused  by  impending  singularity  of  the  operator  (A  -  AB) 
as  the  bifurcation  point  is  approached,  and  is  manifested  by  ill-conditioned 
equations  leading  to  unreliable  values  for  the  changes  in  the  prebuckling 
configuration.  A  technique  has  been  devised  which  actually  uses  the  fact 
that  the  operator  (A  -  AB)  is  singular  to  determine  the  changes  in  the  pre¬ 
buckling  configuration  accurately. 

As  Koiter  (1945)  points  out,  the  eigenvector  is  orthogonal  to 
changes  in  the  prebuckling  configuration  at  the  bifurcation  point.  A  side 
condition  is  thus  available  in  the  form 

XT  C  Y  =  0  (2,7) 


or  in  the  form  of  an  equivalent  integral  expression  when  X  and  Y  are  con¬ 
tinuous  quantities.  The  X  and  Y  refer  to  the  eigenvector  and  incremental 
change  of  the  prebuckling  configuration,  respectively.  The  quantity  C  is 
a  suitable  self-adjoint  positive-definite  operator.  This  device  is  employed 
only  for  the  determination  of  accurate  changes  in  the  prebuckling  configura¬ 


tion  near  the  bifurcation  point. 


The  actual  choice  of  C  rs  indicated  for 


a  particular  example  in  Chapter  4. 

The  addition  of  Eq.  (2.7)  to  the  system  of  equations  to  be  solved 
for  the  incremental  changes  in  the  prebuckling  configuration  means  there 
are  now  more  equations  than  unknowns.  Actually  the  equations  are  not  all 
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independent:  at  the  bifurcation  point.  It  appears  to  be  easiest,  from  a 
computational  standpoint,  to  derive  an  independent  set  of  equations  by 
pre-multiplying  the  equations  by  the  transpose  of  the  coefficient  matrix. 
This  is  equivalent  to  the  so-called  least-squares  technique.  Indeed,  away 
from  the  bifurcation  point,  a  least-squares  interpretation  of  the  com¬ 
putation  is  appropriate  because  the  equations  are  independent.  Appending 
the  side  condition  to  the  original  equations  results  in 

D  y  =  b  (2.8) 

where  D  has  one  more  row  than  column.  The  least  squares  solution  of 
Eq.  (2.8)  yields 

DTD  y  =  DTb  (2.9) 

1' 

For  the  algebraic  eigenvalue  problem,  the  matrix  D  D  may  be  shown  to  be 
nonsingular  (see  Appendix  B) . 

2.5.  The  Initial  Eigenvector 

The  method  of  generating  the  initial  eigenvector  is  most  easily 
explained  in  the  context  of  a  particular  problem  and  solution  technique. 
However,  in  Section  2.2  of  this  chapter,  a  method  of  generating  the  eigen¬ 
vector  for  the  algebraic  eigenvalue  problem  is  outlined  for  the  special 
case  of  A,  B  and  A  corresponding  to  the  onset  of  buckling.  An  approximate 
eigenvector  may  be  generated  in  the  same  way  even  though  A,  B  and  X  do 
not  correspond  to  buckling.  It  has  been  found  that  some  care  must  be  taken 
in  the  process  of  finding  the  approximate  eigenvector.  This  matter  will 
be  discussed  in  detail  in  Chapter  4. 
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2.6.  Observations  and  Comments 

Although  the  technique  is  examined  for  the  cases  when  A  and  B 
depend  on  the  eigervaiue  A ,  it  should  be  evident  that  several  types  of 
less  complicated  eigenvalue  problems  are  encompassed  by  this  general 
theory.  For  instance,  buckling  l^ads  of  Euler  struts  and  the  modes  of 
small-amplitude  free  vibration  of  elastic  systems  are  examples  where  A  and 
B  do  not  depend  on  the  eigenvalue.  In  fact,  the  technique  x.?as  first  tested 
on  these  simpler  problems. 

By  restricting  A  and  B  tc  be  self-adjoint  and  positive-definite, 
it  is  possible  to  place  some  aspects  of  the  proposed  method  on  a  firm 
theoretics;  basis  (see  Appendices  A  and  R) .  In  addition,  physical  argu¬ 
ments  and  experience  in  solving  a  number  of  problems  provide  considerable 
evidence  for  the  wide  applicability  of  the  method. 

A  paper  by  Rail  (1961)  proposed  an  iterative  procedure  for  finding 
eigenvalues  and  eigenvectors  of  a  discrete  system.  There  is  a  formal 
relation  between  Rail's  method  and  the  present  one,  but  in  Rail's  method 
the  eigenvalue  is  not  treated  as  an  unknoxvn  the  same  basis  as  the  components 
of  the  eigenvector.  Further,  in  Rail's  method  there  is  no  freedom  in  the 
choice  of  a  ''side  condition"  and,  in  fact,  an  unfortunate  choice  of  co¬ 
ordinates  can  lead  to  failure  of  the  procedure. 
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3.  THE  PREBUCKLING  CONFIGURATION 


3.1.  Introduction 

In  Chapter  2,  a  general  technique  is  presented  for  the  simul¬ 
taneous  Improvement  of  an  approximate  bifurcation  point  and  eigenvector. 
There  the  technique  is  presented  generally  and,  therefore,  somewhat  ab¬ 
stractly.  In  Chapters  3  and  4  the  solution  process  for  the  buckling  of  a 
rod-type  member  is  presented  in  some  detail  as  an  example  of  the  use  of 
the  general  technique  of  Chapter  2.  The  nature  of  the  technique  requires 
a  method  of  determining  an  equilibrium  configuration  corresponding  to  a 
given  load  level  which  in  general  is  given  by  the  solution  of  a  system  of 
nonlinear  equations.  The  procedure  for  solution  of  the  nonlinear  equations 
at  some  distance  from  a  bifurcation  point  is  presented  in  this  chapter. 

3.2.  Problem  Description 

For  a  detailed  analysis  of  the  arch  problem,  the  equations  ex¬ 
pressing  the  three-dimensional  behavior  of  a  rod-type  member  will  be  pre¬ 
sented  and  their  method  of  solution  described.  Since  the  boundary  con¬ 
ditions  and  loading  are  pertinent  to  the  analysis,  a  specific  choice  must 
be  nisds*  Hsrs  ths  msnibsr  will  bs  sssuinsd  to  bs  clsnipsd  st  tho  boundsriss 
and  loaded  with  a  concentrated  load  (see  Fig.  3(b)  and  Fig.  4). 

As  mentioned  in  Chapter  1,  the  equilibrium,  geometric  and  moment- 
curvature  relationships  are  those  presented  by  Love  (1927).  Love  also 
indicates  how  these  equations  must  be  modified  in  order  to  include  the  ef¬ 
fects  of  extension  of  the  member  centerline  In  this  study,  extension  of 
the  centerline  is  neglected  for  the  full  three-dimensional  problems,  al¬ 
though  results  will  be  presented  in  Chapter  5  for  some  two-dimensional 
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problems  where  extension  of  the  centerline  is  included.  The  effects  of 
restraint  of  warping  of  the  member  cross  section  are.  not  included  in  the 
discussion  of  this  chapter,  but  results  are  presented  in  Chapter  5  for 
lateral  buckling  of  an  initially  straight  I-beam  under  a  dead  load  where 
warping  restraint  is  considered.  Timoshenko  and  Gere  (1961)  and  Vlasov 
(1959)  indicate  the  formulation  of  the  proper  equations  relating  the 
twist  of  the  member  to  the  torsional  moment  when  restraint  of  warping  is 
considered. 

3.3.  Basic  Equations  for  the  Behavior  of  an  Initially  Curved  Member 

3.3.1.  Preliminaries 

Figure  1  shows  the  member  and  global  coordinate  system.  Two  of 
the  member  axes  are  taken  as  the  principal  axes  of  the  section  and  the 
third  is  directed  along  the  tangent  to  the  centerline  of  the  member.  The 
member  and  global  coordinate  systems  are  related  by  the  following  matrix 
transformation. 


where  the  £. ,  m . ,  and  n.'s  are  direction  cosines, 
i  i  r 
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3.3.2.  Equilibrium  Equations 

The  equations  of  equilibrium,  as  presented  by  Love  (1927)  may  be 

written  as 


K.  N,  =  0 
J  k 


dM 

ds~  "  eijk  Kj  \  “  e3ik  Hk  =  ° 


(3.2) 


The  summation  convention  will  be  used  throughout,  unless  the  contrary  is 
specifically  stated.  Also,  the  subscripts  i,  j,  k  will  always  take  on  the 
values  1,  2,  3.  The  quantities  N^,  are  internal  forces,  internal 

moments  ar  curvat  re  vectors,  respectively,  in  the  local  coordinate  system. 
The  quantity  e ...  is  the  alternating  tensor  and  s  is  the  arc  length. 

XJK 

3.3.2.  Geometric  Equations 

Although  there  are  only  three  independent  direction  cosines,  it 
is  convenient  to  ignore  this  fact  temporarily  and  to  present  the  entire  set 
of  geometric  equations.  The  nine  equations,  relating  direction  cosines  to 
curvatures  are 

d£ 

dT  '  eijk  =  0 


dm^ 

ds 


Eljk  =  0 


(3.3) 


dn 


i 


ds 


cijk  “A ' 0 
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3.3.4.  Displacement  Equations 

The  equilibrium  and  geometric  equations  do  not  involve  displace¬ 
ments  explicitly.  However,  the  equations  expressing  satisfaction  of  the 
boundary  conditions  do,  in  general,  involve  displacements.  The  displace¬ 
ment  quantities  required  are  derivable  from  the  direction  cosines  by  a 
simple  quadrature. 

15 

X1(s)  =  f  13(C)  d£ 

J  0 

s 

m3(0  d£  (3.4) 

0 

s 

X3(s)  =  j  n3<0  d5 
0 

where  £  is  a  dummy  variable  and  the  X^(s)  are  the  global  coordinates  of 
the  centerline  of  the  member  as  functions  of  the  arc  length,  s. 

3.3.5.  Moment-Curvature  Relations 

The  effects  of  restraint  of  warping  are  not  considered  in  the 
behavior  of  the  arches  studied  here.  Thus  the  torsional  behavior  is 
entirely  of  the  St.-Venant  type.  The  torque  is  given  by  the  product  of  the 
change  of  the  rate  of  twist,  K3  -  K3q,  and  the  St.-Venant  torsional 
rigidity,  GJ.  For  consistency  of  notation,  GJ  is  taken  equal  to  EI^. 
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Thus  the  moment-curvature  relations  become 

=  EI^  (K^  -  K^,q),  (no  summation,  i  =  1,  2,  3)  (3.5) 

where  EI^  refers  to  the  various  rigidities  and  is  the  curvature  vector 
in  the  unloaded  state. 

3.3.6.  Conditions  at  a  Concentrated  Load 

The  global  representation  of  the  concentrated  load  is  taken  as 

P  =  P  I2  (3.6) 

where  I 2  is  a  unit  vector  in  the  global  X^  direction  and  P  is  the  magnitude 
of  the  force,  which  is  assumed  to  be  applied  at  the  centerline  of  the 
member . 

Consideration  of  equilibrium  of  an  element  of  arch  containing 
the  concentrated  force  yields  the  following  "jump  conditions"  relating  the 
internal  force  resultants  on  either  side  of  the  load, 

+  P  m  -  =  0  (3.7) 

The  superscripts  +,  refer  to  points  to  the  right  and  left  of  the  load, 
positive  being  in  the  direction  of  increasing  arc  length. 

3.3.7.  Boundary  Conditions  for  a  Clamped  Arch 

For  a  clamped  arch,  the  boundary  co:  dltions  specify  that  both 
the  direction  cosines  at  the  supports  and  the  global  coordinates  of  the 
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supports  remain  unchanged.  The  boundary  conditions  for  an  initially  planar 
clamped  arch  are 


/  =  £ 

2  20 


ml  "  m10 


n3  =  n30 


Xi  =  *10 


(at  s  =  0,  s  =  sf) 


(3.8) 


/here  the  second  subscript  0  indicates  the  original  configuration  and  s^  is 
the  arc  coordinate  of  the  far  boundary. 

3.3.8.  Complementary  Loading  Parameter 

It  has  been  noted  previously  by  Bueckner,  Johnson  and  Moore  (1965) 
and  Leicester  (1068)  that  a  numerical  analysis  of  snap-through  buckling  of 
shallow  spherical  shells  can  encounter  difficulties  associated  with  the  in¬ 
cremental  loading  process.  A  similar  difficulty  occurs  in  arches.  This 
difficulty  stems  from  the  fact  that  so-called  limit  points  (see  Fig.  2)  may 
exist  m  the  force-deflection  curve.  If,  near  point  A  an  increment  of  force 
is  chosen  such  that  the  total  force  is  greater  than  P^,  obviously  there 
is  no  solution.  This  Is  a  very  real  possibility  since  in  general  the  maximum 
value  is  not  known  in  advance.  Near  point  A,  the  difficulty  may  be  over¬ 
come  by  incrementing  the  deflection  instead  of  the  force.  A  similar  situa¬ 
tion  occurs  near  point  B  except  that  the  force  quantity  should  be  incremented 
instead  of  the  deflection.  In  the  vicinity  of  the  limit  points,  convergence 
of  the  Newton-Raphson  or  successive  approximation  procedures  will  be  slow 
or  fail  entirely  if  a  poor  choice  of  loading  parameter  is  made.  For  this 
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reason  it  is  advantageous  to  be  able  to  select  either  force  or  deflection 
as  the  independent  variable  in  the  loading  process. 

In  order  to  demonstrate  how  a  loading  parameter  other  than  the 
concentrated  force  itself  is  used  in  the  solution  process,  a  complementary 
loading  parameter  corresponding  to  the  deflection  under  the  concentrated 
force  and  in  the  direction  of  the  force  will  be  used  here.  The  expres¬ 
sion  for  this  component  of  the  deflection  under  the  load  is 


s 

d2  =  f  “  m30(5)}  d5  (3'9) 

J  0 

where  the  upper  limit  of  integration,  s^,  refers  to  the  arc-length  coor¬ 
dinate  of  the  point  of  application  of  the  force. 

3.4.  Solution  of  Nonlinear  Equations 

3.4.1.  General  Discussion 

There  are  several  techniques  available  for  solving  cwo-point 
boundary  value  problems  described  by  nonlinear  ordinary  differential  equa¬ 
tions.  The  character  of  t"  2  particular  set  of  equations  may  limit  the 
effectiveness  of  some  of  these  techniques . 

One  particular  technique  called  the  "shooting  method"  has  been 
used  by  Huddleston  (1968)  to  solve  the  nonlinear  equations  which  describe 
the  large  deflections  of  an  arch  under  a  concentrated  load.  The  boundary 
value  problem  is  converted  to  an  initial  value  problem  and  the  nonlinear 
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equations  integrated  numerically .  Since  some  of  the  initial  values  are 
unknown,  these  are  adjusted  until  the  far  boundary  conditions  are  satisfied. 
Generally  a  few  iterations  are  required  to  satisfy  the  boundary  conditions. 
This  technique  will  encounter  numerical  difficulties  when  the  solution  of 
the  nonlinear  ordinary  differential  equation  can  exhibit  a  boundary  layer 
or  edge  effect.  In  this  case,  the  initial  value  solutions  will  grow  rapidly 
as  they  are  propagated  to  the  far  boundary.  Since  computers  carry  a  finite 
number  of  digits  in  numerical  computations,  the  quantities  required  for  the 
equations  which  express  satisfaction  of  the  far  boundary  conditions  may  have 
literally  no  significance  because  of  round-off  during  the  numerical  inte¬ 
gration  process.  In  fact,  this  phenomenon  can  occur  even  though  the  initial 
values  are  quite  close  to  the  correct  ones. 

Another  technique  has  been  developed  by  Berezin  and  Zhidkov  (I960) 
and  by  Jordan  and  Shellay  0.966'  for  solving  just  the  ‘ype  of  problem  where 
"growing"  solutions  are  present.  This  technique  does  not  require  iteration 
but  a  transformation  of  the  equations  to  a  new  set  of  /ariables  is  necessary 
before  the  solution  may  proceed.  As  with  the  "shooting  method",  the  trans¬ 
formed  set  of  equations  is  in:cgrated  numerically  since  they  are  in  general 
nonlinear.  Jordan  and  Shelley  indicate  that  if  the  original  problem  does 
not  have  a  bouidary  or  eage  effect,  the  transformed  solution  may.  In  this 
case,  the  transferee i  problem  would  encounur.  numerical  difficulties.  It 
turns  out  that  even  if  there  is  a  boundary  effect,  it  is  possible  that  the 

■k 

method  will  fail. 

The  technique  used  in  this  study  does  not  depend  on  the  character 
of  the  nonlinear  equations.  That  is,  the  presence  of  a  boundary  or  edge 


* 


This  observation  ±b  due  to  Professor  M.  S.  Zarghamee. 
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effect  does  not  present  any  serious  obstacles.  The  Newton-Raphson  technique 
is  used  to  solve  the  nonlinear  equations  and  thus  only  linearized  equations 
are  integrated.  When  growing  solutions  are  present  in  the  integration  of 
the  linearized  equations,  the  suppression  technique  used  by  Zarghamee  and 
Robinson  (1967)  and  Goldberg,  Setlur  and  Alspaugh  (1965)  is  implemented  to 
avoid  the  loss  of  significant  figures  due  to  round-off. 

3.4.2.  The  Newton-Raphson  Procedure 

The  nonlinear  equations  of  this  study  are  solved  using  the  Newton- 
Raphson  procedure.  In  the  use  of  this  procedure,  the  loading  is  applied  to 
the  structure  in  increments  (not  necessarily  small)  by  the  following  com¬ 
putational  process.  The  reason  for  applying  the  loading  parameter  in  steps 
will  become  apparent  as  the  discussion  proceeds. 

Assume  that  at  some  jtage  in  the  loading  process  a  solution  E-^ 
of  the  nonlinear  equations  is  known  which  corresponds  to  a  loading  level  L^. 
An  increment  of  load  AL^  is  now  applied,  i’he  Newton-Raphson  procedure  is 
used  to  find  a  new  equilibrium  configuration  corresponding  to  the  total 
loading  parameter  given  by  L^.  +  AL^.  The  equations  specifying  the  linear 
response  of  the  configuration  must  then  be  derived  by  linearizing  the 
equations  about  this  configuration.  '’’he  linear  incremental  solution  6Ej. 
corresponding  to  an  increment  of  loading  AL^.  is  added  to  the  existing  con¬ 
figuration  E  to  produce  a  new  configuration  E_  , .  In  general  the  con- 
X  i“rX 

figuration  will  not  satisfy  the  nonlinear  equations  since  a  linear 
approximation  was  used  to  compute  <5  E^.  Thus  there  are  residuals  in 
these  nonlinear  equations. 
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The  next  step  is  to  remove  the  residuals,  without  a  further  increase 
in  the  loading  parameter.  The  equations  are  again  linearized,  this  time 
about  the  new  configuration  The  linear  response  o  E^.  at  c^is  c''n" 

figuration  is  calculated.  The  "loading"  in  this  computation  consists  of 
the  negatives  of  the  residuals  in  the  corresponding  nonlinear  equations.  A 
new  configuration  E  j  squal  to  ^  is  thus  derived.  At  this  point, 

the  configuration  is  substituted  into  the  nonlinear  equations  and  the 

resulting  residuals  are  again  examined.  If  the  residuals  are  t,mall  enough, 
a  new  equilibrium  configuration  has  been  found  and  another  increment  of  the 
loading  parameter  may  be  applied.  If  the  residuals  are  not  satisf 'otcry , 
this  process  cf  removing  residuals,  for  a  constant  value  of  loading  para¬ 
meter,  is  repeated  until  a  new  equilibrium  configuration  is  obtained. 

It  is  evident  from  the  above  discussion  that  it  is  necessary  to 
linearize  the  nonlinear  equations  of  Sections  3.3.2.  -  3.3.8.  about  a 
general  reference  configuration  in  order  to  use  the  Newton-Raphson  procedure. 
These  linearized  equations  are  presented  in  the  next  section. 


3.4.3.  Linearization  of  the  Prebuckline  Configuration 


In  order  to  avoid  the  cumbersome  notation  of  Chapter  2  in  expres¬ 
sing  the  linearized  equations  of  the  arch  problem,  the  superscript  j  used 

tH 

in  Chapter  2  to  denote  the  j  configuration  will  ^e  dropped  and  instead 
the  current  configuration  will  instead  be  denoted  simply  by  the  quantities, 
Ki*  ^i*  ^i»  mi?  nis  etc*  a  superscript.  Since  the  equations 

specifying  the  prebuckling  configuration  are  of  first-order,  the  lineariza¬ 


tion  process  is  particularly  straightforward  and  leads  to  the  following 
equations . 


i 
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Linearized  Equilibrium  Equations: 


6 (dN±) 
as 


EiJit  ({Kj  \  +  5\>  “  0 


SUM,) 

S ~  -  ei]k  (SKj  %  T  KJ  SMk>  -  e3ik  6Kk  "  0 


(3.10) 


Linearized  Geometric  Equitions: 


6(d££) 

ds 


e..,  (6Kjn  +  K  .St) 
ijk  J  1c  j  K 


0 


6  (dm . ) 

“di”  -  eijk  (6Kj  *k  +  Kj  5\)  =  ° 


(3.11) 


S(dn.) 

ST  -  eijk  (6Kj  °k  +  Kj  5V  =  0 


Linearized  Displacement  Equations: 


SX1  =  J  (5)  dC 


/ 


CX2  =  j  Sm3  (O  d* 


(3.12) 
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C3  =  f  6i>3  (O 

0 


Linearized  Moment-Curvature  Relations: 


5M^  =  (21)^  6K^  ,  (no  summation) 


Linearized  Condition  at  the  Concentrated  Load: 


6N,)  '  +  P  6  m  +  6  P  mi  =  0 


Linearized  Boundary  Conditions: 


6l2  -  0 


<5m^  =  0 


6n3  =  0 


6X±  =  0 


(at  s  =  0,  s  =  sf) 


Linearized  Complementary  Loading  Parameter: 


■/ 


$n3  (5)  d5 


(3.12) 


(3.13) 


(3.14) 


(3.15) 


(3.16) 


The  6N. ,  dM^>  dlC^,  dm^,  dn^,  dd3,  etc.,  are  the  linearized 

quantities  where  tne  d  is  used  to  denote  a  linear  increment.  In  general, 
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Eqs.  (3.10),  (3.11),  (3.12),  (3.14),  (3.15),  (3.16)  when  they  are  applied, 
will  have  on  their  right  hand  sides  not  zeros  but  the  negatives  of  the 
residuals  computed  from  their  corresponding  nonlinear  equations  as  explained 
in  Section  3.4.2. 

3.5.  Typical  Incremental  Loading  Cycle 

The  typical  incremental  loading  cycle  of  this  study  may  be  sum- 
maried  as  follows  using  the  notation  which  has  been  introduced: 

(1)  Assume  that  an  equilibrium  configuration  corresponding 

to  the  quantities  M^,  ,  K.^,  Z^,  m^,  n^,  d^,  etc.  is 

known; 

(2)  Apply  an  increment  &  d2  of  the  loading  parameter  by  use 
of  the  linearized  equations  (Eqs.  (3.10)  -  (3.16))  to 
obtain  6N_^,  oM^,  6K^,  &Z^,  <$m^,  6n^,  Ac^,  etc.; 

(3)  Add  the  incremental  quantities  6N_^,  <5K^,  SZ^,  6m^, 

<$ni’  Ad2>  etc.  determined  in  the  previous  step  to  the 
previous  values  of  N^,  K^,  Z^,  nu,  n^,  d2»  etc.  to 
obtain  a  new  set  of  N^,  M^,  K_^,  Z^,  ,  n^.,  d2>  ecc.; 

(4)  Compute  the  residuals  in  Eqs.  (3.2)  -  (3.9)  using 

the  new  N^,  Z^,  in, ,  n^,  d2>  etc.  of  step  (3); 

(5)  Check  the  residuals  to  see  if  they  are  acceptable.  If 
so,  the  process  stops,  a  new  equilibrium  configuration 
having  been  determined.  If  the  residuals  are  not 
acceptable,  go  on  to  step  (6),  Note  there  are,  in 
general,  residuals  in  the  jump  condition  Eq.  (3.7)  and 
in  the  complementary  loading  parameter  expression 
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Eq.  (3.9)  as  well  as  the  differential  equations; 

(6)  Remove  the  residuals  obtained  in  step  (4)  by  computing 
the  linear  effect  on  the  new  configuration  (deter¬ 
mined  in  step  (3))  of  the  negatives  of  the  residuals 
determined  in  step  (4).  Go  back  to  step  (3). 

Although  the  same  equations  are  used  in  steps  (2)  and  (6) , 

(except  for  the  right  hand  sides)  the  increase  in  the  loading  parameter  d^ 
is  carried  out  only  once.  Note  that  the  N^,  M^,  K^,  nu,  n^,  d^,  etc. 
are  always  the  latest  quantities. 

3.6.  Details  of  the  Solution  of  the  Linearized  Differential  Equations 

The  discussion  of  a  typical  incremental  loading  cycle,  Section  3.5, 
was  based  on  the  assumption  that  a  solution  to  the  two-point  boundary  value 
problem  given  by  the  linearized  differential  equations,  boundary  conditions, 
jump  condition  and  incremental  loading  parameter,  could  be  found.  In  this 
study,  the  modified  two-point  boundary  value  problem  defined  by  the  linearized 
differential  equations,  the  boundary  conditions,  jump  condition  and  the  in¬ 
cremental  loading  equation  is  converted  to  an  initial  value  problem.  The 
initial  value  technique  has  been  used  by  Kalnins  (1964),  Gcldbci.6,  Setlur, 
and  Alspaugh  (1965),  and  Zarghamee  and  Robinson  (1967)  to  solve  boundary 
value  problems  described  by  ordinary  differential  equations.  Since  the 
method  uses  one  boundary  as  the  origin  of  the  linearized  initial  value 
problem,  the  so-called  initial  values  are  selected  so  as  to  satisfy  the 
boundary  conditions  at  the  origin  automatically.  As  the  method  is  used  here, 
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a  set  of  inderandent  initial  value  solutions  (see  Table  1)  is  propagated 
from  the  origin  to  the  far  boundary  where  a  linear  combination  of  these 
soli  _’'ons  is  formed  to  satisfy  the  jin.carized  boundary  conditions  and  the 
condition  on  the  incremental  loading  parameter  Eq.  (3.16). 

The  increments  in  the  boundary  displace  at  the  far  end  and 

in  the  loading  parameter  are  expressed  as  integrals  of  the  quantities  oc¬ 
curring  in  the  linearized  differential  equations.  This  means  that  the 
equations  (incremental  boundary  conditions  and  incremental  loading  para¬ 
meter)  for  determining  the  proper  linear  combination  of  solutions  require 
that  a  quadrature  of  the  quantities  in  the  individual  initial  value  solu¬ 
tions  be  carried  out.  This  has  been  done  numerically  using  Simpson's  rule. 
The  condition  on  the  incremental  loading  parameter  is  treated  the  same  as 
an  additional  boundary  condition  when  forming  the  linear  combinations  neces¬ 
sary  to  solve  for  the  correct  initial  values. 

The  individual  initial  value  solutions  are  found  by  n'merical  in¬ 
tegration  using  a  trapezoidal  integration  formula  as  part  of  a  predictor- 
corrector  process.  Tl.e  numerical  integration  process  has  been  presented  by 
Crandall  (1956) .  The  character  of  these  equations  is  such  that  rapidly 
growing  solutions  are  not  present  in  the  numerical  integration  process.  For 
this  reason,  the  so-called  suppression  technique  (see  Section  3.4.1.)  is 
not  necessary.  In  Chapter  5  of  this  study,  an  example  problem  of  the 
lateral  buckling  of  an  I-beam  with  warping  rigidity  is  solved  which  requires 
suppression  during  the  integration  process. 

Table  1  shows  the  initial  values  for  each  solution.  The  residual 
terms  in  the  particular  solution  occur  in  what  has  been  called  step  4  of  the 
incremental  loading  procedure  as  given  in  Section  3.5. 
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3.7.  Direction  Cosine  Correction 

Since  the  direction  cosines  are  treated  as  independent  quantities 
during  the  numerical  integration  of  the  linearized  differential  equations, 
it  is  possible  that  "drift"  of  the  direction  cosines  will  take  place  so  that 
they  will  no  longer  form  an  orthonormal  set.  The  computational  process 
guarantees  that  the  squares  and  scalar  products  of  the  new  local  coordinate 
basis  vectors  are  constant  across  the  arch.  However  there  is  no  mechanism 
in  the  straightforward  procedure  to  control  drift  in  these  constants,  which 
should,  of  course,  be  either  one  or  zero.  A  technique,  outlined  in  Appendix 
C,  has  been  developed  to  ensure  orthonormality. 

3.8.  Other  Boundary  Conditions 

An  arch  which  is  simply  supported  in  the  plane  presents  no  added 
complications.  The  geometric  boundary  condition,  n^  =  n,g  is  replaced  by 
the  moment  condition  =  0.  See  Figs.  3(a)  and  3(b). 

Other  types  of  Loundary  conditions  may  require  considerable  care 
in  their  formulation.  If  it  is  desirable  to  allow  more  than  one  free 
rotation  at  a  support,  it  is  useful  to  have  in  mind  a  physical  model  (say 
a  Hooke's  joint)  of  the  support  in  order  to  avoid  the  possibility  of  intro¬ 
ducing  a  nonconservative  force  system  at  the  support.  This  difficulty  has 
been  explained  in  detail  by  Ziegler  (1956) . 
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4.  DETERMINATION  OF  POINTS  OF  BIFURCATION  IN  THE  CASE  OF 
NONLINEAR  PREBUCKLING  BEHAVIOR 

4.1.  Introduction 

As  mentioned  in  Chapter  2,  a  study  of  postbuckliig  brhavior  re¬ 
quires  the  location  of  the  bifurcation  point.  This  chapter  deals  with  a 
specif' c  application  of  the  general  technique  of  Chapter  2  for  improving 
an  approximation  to  a  bifurcation  point  and  the  corresponding  approximate 
eigenvector.  For  the  specific  arch  problem,  a  prebuckling  coni iguration 
determined  by  the  method  of  Chapter  3  is  used  as  an  approximation  to  the 
bifurcation  point  in  the  process  described  in  Chapter  2.  The  method  for 
generating  the  corresponding  approximate  eigenvector  will  be  given  in  detail 
later  in  this  chapter.  Since  this  technique  requires  not  only  a  knowledge 
of  the  local  behavior  of  the  prebuckling  configuration  (the  Y  of  Sec.  2.2.) 
but  also  the  eigenvector  "branching"  from  a  prebuckl:ng  curve  (the  X  of 
Sec.  2.2.).  two  different  incremental  quantities  must  be  studied  at  the 
same  time.  It  is  not  difficulc  to  adapt  the  linearized  equations  of  Chapter 
3  for  this  purpose  with  a  suitabi’  change  of  notation.  The  new  linearized 
equations  will  be  solved  for  the  quantities  corresponding  to  the  eigenvector, 
which  is  "along"  the  initial  segment  of  a  new  branch.  These  linearized 
equations  will  be  referred  to  as  "branch  equations". 

4.2.  Branch  Equations 

The  following  equations  are  the  linearized  equations  of  Chapter  3 
with  the  4  replaced  by  an  asterisk.  As  the  discussion  proceeds,  it  will  be 
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obvious  that  a  new  notation  is  necessary  for  clarity.  These  equations 
play  the  role  of  Eq.  (2.1)  of  Chapter  2. 


Branch  Equilibrium  Equations: 


dNj 

ds 


£«k  (iS  +  Ej  <>  '  0 


* 

dM. 

T_ 

ds 


eijk  Mk  +  Kj  V  "  e3ik  \  = 


=  0 


(4.1) 


Branch  Geometric  Equations: 


ds 


eijk  (Kj  *k  +  Kj  V 


=  0 


k  k 

ds~  '  eijk  (Kj  “k  +  Kj  \) 


(4.2) 


dn. 


as 


eijk  (Kj  °k  +  Kj  V  ° 


Branch 'Displacement  Equations: 


X*(s) 


/ 


l3  (C)  d5 


X*(s) 


J 


m^  (O  di, 


(4.3) 
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X*(S) 


I 


n3  (?)  d? 


(4.3) 


Branch  Moment-Curvature  Relations : 


=  (EI)^  ,  (no  summation) 


(4.4) 


Branch  Condition  at  the  Concentrated  Load: 


*(+)  *(_)  * 

N. V  1  -  N/  J  +  P  m.  =  0 

x  i  x 


(4.5) 


Branch  Boundary  Conditions: 


C  -  0 


m,  =  0 


-  0 


Xi  "  0 


(at  s  ■=  0,  s  =  sf) 


(4.6) 


If  the  prebuckling  configuration  given  by  the  quantities  N^, 

K^,  m^,  n^,  etc.  is  the  one  corresponding  to  bifurcation,  the  eigenvector 

may  be  generated  from  Equations  4.1  -  4.6  in  a  straightforward  manner.  In 
general  this  fortuitous  circumstance  will  not  previil  and  the  prebuckiing 
configuration  must  be  adjusted  in  order  to  reach  the  bifurcation  point.  The 
crux  of  the  problem  then  ic  to  adjust  the  prebuckling  configuration  so  that 
a  better  approximation  to  the  bifurcation  point  is  obtained.  The  general 
technique  developed  in  Chapter  2  is  used  for  this  purpose. 


37 


As  ,ume  that  an  approximate  prebuckling  configuration  found  by 
the  method  of  Chapter  3  and  an  approximate  eigenvector  are  substituted  into 
Eqs.  (4,1)  -  (4.6).  There  is,  in  general,  a  residual  ir.  these  equations. 
The  modification  of  the  Newton-Raphson  technique  introduced  in  Chapter  2 
is  used  to  remove  the  residuals.  Here  it  is  necessary  to  linearize  the  so- 
called  branch  equations  with  respect  to  the  prebuckling  (unstarred)  quan¬ 
tities  (N^,  M^,  K^,  4L,  nu,  n^,  etc.)  and  the  current  approximate  eigen- 

******* 
vector  (N^,  nu,  n^,  etc.). 

4.3.  Linearized  Branch  Equations 

As  noted  in  Chapter  2,  two  types  of  incremental  quantities  appear 
in  the  linearized  branch  equations;  those  corresponding  to  changes  of  the 

prebuckling  configuration  (6N^,  61^,  6nu,  6n^,  etc.)  and  those 

*  *  *  .*  *  .  * 
corresponding  to  changes  in  the  eigenvector  (6N^,  6.L,  dK^,  6-t^,  6nu,  6n^, 

etc.).  The  linearized  branch  equations  are  understood  to  be  valid  about 

a  "hyper-configuration"  consisting  of  the  current  prebuckling  configuration 

and  the  approximate  eigenvector.  Also,  in  general,  Eqs.  4.7,  4.8,  4.9, 

4.11,  and  4.12  will  have  non-zero  right  hand  sides  equa]  to  the  residuals 

computed  from  the  corresponding  nonlinear  branch  equation.  The  linearized 

branch  equations,  as  given  below,  play  the  role  of  Eq.  (2.3). 

Linearized  Branch  Equilibrium  Equations: 

* 

<$dN.  ^  A  A  * 

1  -  (SKj  \  *  Kj  4\  +  5  Kj  "k  +  Kj 5  V  ■ 


ds 


0  (4.7) 
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As  »ume  that  an  approximate  prebuckling  configuration  found  by 
the  method  of  Chapter  3  and  an  approximate  eigenvector  are  substituted  into 
Eqs.  (4,1)  -  (4.6).  There  is,  in  general,  a  residual  ir.  these  equations. 

The  modification  of  the  Newton-Raphson  technique  introduced  in  Chapter  2 
is  used  to  remove  the  residuals.  Here  it  is  necessary  tc  linearize  the  so- 
called  branch  equations  with  respect  to  the  prebuckling  (unstarred)  quan¬ 
tities  (N^,  M^,  K^,  m^,  n^,  etc.)  and  the  current  approximate  eigen- 

^  ^  ^  ^  A  ^ 

vector  (N. ,  M . ,  K.,  t.  ,  m. ,  n.,  etc.), 
i  i  x  x  l  l 

4.3.  Linearized  Branch  Equations 

As  noted  in  Chapter  2,  two  types  of  incremental  quantities  appear 

in  the  linearized  branch  equations;  those  corresponding  to  changes  of  the 

prebuckling  configuration  (6N^,  61^,  6Z^,  6nu,  6n^,  etc.)  and  those 

^  ^ 

corresponding  to  changes  in  the  eigenvector  (dN^,  6.L,  6K^,  6Z^,  6nu,  6m, 
etc.).  The  linearized  branch  equations  are  understood  to  be  valid  about 
a  "hyper-configuration"  consisting  of  the  current  prebuckling  configuration 
and  the  approximate  eigenvector.  Also,  in  general,  Eqs.  4.7,  4.8,  4.9, 

4.11,  and  4.12  will  have  non-zero  right  hand  sides  equal  to  the  residuals 
computed  from  the  corresponding  nonlinear  branch  equation.  The  linearized 
branch  equations,  as  given  below,  play  the  role  of  Eq.  (2.3). 

Linearized  Branch  Equilibrium  Equations: 

6dN*  *  *  *  * 

IT  -  £«k  <{Kj  "k  -  Kj  “k  +  6  \  +  Kj  «  V  -  0  W.7) 


0 
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Linearized  Branch  Moment-Curvature  Relations: 

*  (EI)^  ,  (no  summation)  (A. 10) 


Linearized  Branch  Condition  at  the  Concentrated  Load: 


•%(+)  *(-)  *  * 

6N.  -  6N.  '  +  r5m.  +  6Pm,  =  0 

i  i  i  i 

Linearized  Branch  Boundary  Conditions: 

-  0 

6m*  -  0 

(a"  s  =  0,  s  =  sf ) 

* 

5p  =0 

J 

-k 

6Xj,  =  0 


(A. 11) 


(A. 12) 


Since  the  linearized  branch  equations  contain  incremental  terms 
associated  with  changes  of  the  prebuckling  configuration  (the  unstarred 
quantities)  a  preliminary  computation  is  necessary  before  the  actual  solution 
can  proceed.  This  computation  involves  the  determination  of  the  linearized 
response  of  the  prebuckling  configuration  for  6P  =  1;  i.e.,  the  counterpart 
here  of  the  computation  in  Section  2. A.  The  method  for  carrying  out  this 
part  of  the  solution  of  the  linearized  branch  equations  depends  cn  bow 
"close"  the  current  prebuckling  configuration  is  to  the  bifurcation  point. 
Section  A. 5  is  devoted  to  this  topic. 
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1't  is  also  necessary  to  compute  an  initial  approximation  to  the 
eigenvector  before  solving  the  linearized  branch  equations,  as  it  is  the 
interaction  of  the  approximate  eigenvector  with  the  prebuckling  configura¬ 
tion  that  produces  the  residuals  which  "drive''  the  linearized  branch 
equations.  The  computation  of  the  approximate  eigenvector  is  discussed 
in  Section  4.6. 

If  the  approximate  prebuckling  configuration  is  far  enough  from 
the  bifurcation  point  to  permit  use  of  the  standard  Newton-Raphson  technique 
for  the  purpose  of  obtaining  changes  in  the  prrbuckling  configuration,  then 
the  process  of  improving  the  eigenvalue  and  eigenvector  is  straightforward. 
The  linearized  branch  equations  would  form  a  two-point  boundary  value  problem 
except  for  the  fact  that  6P  is  unknown  also.  The  increments  of  the  unstar¬ 
red  quantities  and  6P  are  the  only  unknowns.  The  extra  unknown  6P  is  to 
be  expected  since  the  amplitude  of  the  eigenvector  is  indeterminate.  In 
order  to  solve  the  system  of  linearized  branch  equations,  fi  scalar  side 
condition  is  appended  to  these  equations.  This  side  condition  is  taken  as 


M*  6K*  ds  =  0  (4.13) 

i  i 

0 

This  expression  ensures  that  there  are  not  large  changes  "parallel"  to  the 
eigenvector  when  the  eigenvector  is  close  to  its  true  "direction". 

The  solution  of  these  linearlized  branch  equations  (with  Eq.  (4.13)) 


is  quite  similar  to  the  solution  of  the  linearized  equations  of  Chapter  3. 
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The  scalar  side  condition  introduced  here,  Eq.  C4.13),  plays  the  role  of 
the  complementary  loading  parameter  of  Chapter  3.  The  modified  boundary 
value  problem  described  by  Eqs.  (4.7)  -  (4.12)  and  (4.13),  is  converted  to 
an  initial  value  problem.  As  in  Chapter  3,  a  set  of  initial  value  problems 
is  propagated  from  the  origin  to  the  far  boundary  where  a  linear  combination 
of  these  solutions  is  formed  to  satisfy  the  boundary  conditions  and  the 
scalar  side  condition.  The  procedure  is  similar  enough  to  that  of  Chapter  3 
that,  in  fact,  the  same  numerical  integration  routine  can  be  used  in  both 
parts  of  a  computer  program  to  solve  the  problem.  The  sets  of  initial 
values  given  in  Table  1  carry  over  to  the  solution  process  here  with  the 
understanding  that  the  incremental  branch  quantities  are  now  the  unknowns. 

An  essential  feature  in  the  solution  of  the  linearized  branch  equations  is 
the  presence  of  the  incremental  terms  corresponding  to  changes  of  the  pre¬ 
buckling  configuration.  These  terms  appear  only  in  the  initial  value  solution 
corresponding  to  6P  =  1  (see  Table  1) .  This  should  be  apparent  since  the 
prebuckling  configuration  can  change  only  when  P  changes. 

Once  the  value  of  <5P  is  computed,  the  correct  linearized  change 
in  the  prebuckling  configuration  is  easily  found  by  scaling  the  changes 
caused  by  <5P  =  1,  which  are  found  in  Section  4.5. 

Thus,  both  the  prebuckling  configuration  and  the  eigenvector  are 
modified  simultaneously. 

4.5.  Modifying  the  Prebuckling  Configuration  in  the  Vicinity  of  a 
Bifurcation  Point 

As  indicated  in  Chapters  2  and  3,  there  are  computational  dif¬ 


ficulties  associated  with  computing  the  linearized  response  of  the  prebuckling 


42 


configuration  accurately  in  the  vicinity  of  bifurcation  points.  This 
section  is  devoted  to  a  discussion  of  the  solution  to  this  problem. 

The  changes  in  the  prebuckling  configuration  are  required  to  be 
orthogonal  to  the  eigenvector  (see  Section  2.4).  For  an  inextensional 
centerline,  this  orthogonality  relation  is  conveniently  expressed  as 

Sf 

J  mJ  6K±  ds  =  0  (4.14) 

0 

& 

The  used  in  Eq.  (4.14)  are  the  latest  values  obtained  durir.g  the  process 
of  improving  the  bifurcation  point  and  eigenvector.  This  extra  condition  is 
then  appended  to  the  initial  value  problem  described  in  Chapter  3.  There 
are  now  more  equations  than  unknowns,  but  as  mentioned  in  Chapter  2,  all  of 
these  equations  are  valid  at  the  bifurcation  point.  A  consistent  set  of 
equations  is  derived  using  the  least-squares  technique. 

This  technique  permits  the  accurate  computation  of  changes  in  the 
prebuckling  configuration  near  the  bifurcation  point.  Note,  however,  that 
this  device  is  essential  only  in  the  vicinity  of  the  bifurcation  point.  At 
other  points,  the  standard  Newton-Raphson  technique  outlined  in  Chapter  3  is 
satisfactory  for  modifying  the  prebuekHng  configuration. 

4.6.  Generating  the  Approximate  Eigenvector 

The  process  of  improving  an  eigenvalue  involves  the  solution  of 
a  system  of  non-singular  linear  algebraic  equations.  The  only  difficulty 
is  in  arriving  at  a  suitably  "close"  initial  P  and  eigenvector.  Since  the 
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P  used  is  only  approximate,  there  will  in  general  not  exist  a  solution  of 
the  branch  equations  satisfying  all  the  boundary  conditions.  The  computa¬ 
tional  device  which  has  been  adapted  here  is  to  release  one  of  the  boundary 
conditions.  In  the  first  subsequent  improvement  of  the  P  and  eigenvector, 
it  is  a  straightforward  matter  to  reimpose  the  constraint  which  has  been 
released. 

It  is  obvious  that  there  will,  in  general,  be  more  than  one 
choice  of  constraint  which  can  be  released  for  calculation  of  the  initial 
approximation  of  the  eigenvector.  It  has  been  found  that  by  an  unfortunate 
choice  of  release  of  constraint,  it  is  possible  to  "skip"  the  eigenvalue 
being  sought  and  "jump"  to  a  distant  one.  The  technique  used  to  avoid  this 
problem  is  to  relax  what  appears  to  be  the  "softest"  of  the  constraints. 

For  instance,  in  out-of-plane  buckling  of  ar.  arch,  the  restraint  corresponding 
to  rotation  about  the  tangent  to  the  centerline  at  one  end  of  the  member  is 
relaxed. 

In  general,  it  might  be  necessary  to  run  through  all  choices  of 
constraint  release  at  one  end  to  find  the  one  leading  to  the  smallest  6P 
on  the  first  cycle  of  improvement.  However,  this  extra  computation  is 
actually  not  extensive. 

4.7.  Summary  of  the  Typical  Computational  Cycle 

The  first  part  of  the  cycle  is  really  a  preparatory  stage.  The 
change  in  the  prebuckling  configuration  fcr  <5P  =  1  is  computed  and  the 
approximate  eigenvector  is  generated.  Computational  details  are  explained 
in  Sections  4.5  and  4.6.  At  this  point,  the  current  prebuckling  configuration 
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and  the  approximate  eigenvector  are  substituted  into  the  branch  equations 
and  residuals  are  computed.  These  residuals  are  used  to  "drive"  the 
linearized  branch  equations. 

Because  of  the  way  the  approximate  eigenvector  is  generated, 
during  the  first  iteration  step  the  residuals  do  not  appear  in  the  dif¬ 
ferential  equations  but  only  in  the  boundary  condition  which  was  violated 
when  the  approximate  eigenvector  was  generated.  For  subsequent  iterations, 
there  are,  in  general,  residuals  in  both  the  differential  equations  and 
he  boundary  conditions. 

Eventually,  as  successive  prebuckling  configurations  are  pre¬ 
dicted  and  examined  for  the  presence  of  an  eigenvector,  the  value  of  6P 
and  the  residuals  in  the  branch  equations  computed  during  this  sequence 
'ill  become  acceptably  small.  At  this  point,  the  bifurcation  load  has  been 
reached  and  the  corresponding  eigenvector  generated. 

The  special  process  for  obtaining  changes  in  the  prebuckling 
configuration  when  the  standard  Newton-Raphson  technique  falls  because  of 
poorly  conditioned  equations  was  never  needed  until  the  latest  relative 
change  in  P  was  less  than  U.1C. 

4.8.  Postbuckling  Paths 

Without  referring  to  the  question  of  stability  of  the  paths,  it  is 
a  simple  matter  now  to  proceed  onto  the  branch  given  initially  by  the 
eigenvector.  This  is  done  by  adding  a  multiple  of  the  eigenvector  to  the 
prebuckling  configuration  and  then  dr-ermining  a  new  nonlinear  configura¬ 
tion  using  the  technique  of  Chapter  3. 
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Koiter  (1945)  indicates  that  if  there  is  a  single  branch  from 
the  fundamental  or  prebuckling  path,  stability  of  the  new  path  is  determined 
by  whether  the  load  capacity  increases  or  decreases.  If  the  loaa  in¬ 
creases,  the  new  path  is  stable  and  if  the  load  decreases,  the  new  path  is 
unstable. 

If  there  is  a  multiple  eigenvalue  and  multiple  branches  from  the 
fundamental  branch t  the  stability  considerations  are  more  complicated. 

Koiter  (1945)  has  a  discussion  of  this  more  difficult  problem.  In  Appen¬ 
dix  A  of  this  study,  a  solution  of  the  computational  problem  of  determining 
multiple  branches  is  indicated. 
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5.  NUMERICAL  RESULTS  OF  THE  APPLICATION  OF  THE  THEORY 
TO  ARCHES  AND  BEAMS 


5.1.  General  Remarks 

In  this  chapter,  several  sample  problems  of  the  buckling  of 
arches  are  presented.  In  addition,  a  few  results  are  presented  for  lateral 
buckling  of  a  beam.  These  problems  are  solved  using  the  technique  intro¬ 
duced  in  Chapters  2,  3,  and  4.  The  chief  object  of  these  examples  is  to 
demonstrate  some  of  the  possibilities  of  the  technique.  Compa'"' sons  with 
previous  work  are  made  where  such  work  is  available. 

The  examples  given  in  Sections  5.3.2.  and  5.3.3.  are  planar  arches 
which  may  buckle  only  in  the  plane  of  the  arch  (see  Fig.  (5(b)).  Two  sets 
of  boundary  conditions  and  two  sets  of  rise-to-span  ratios  are  considered. 

In  Section  5.3.4.,  three-dimensional  buckling  of  initially  planar  arches 
is  considered.  That  is,  the  arches  may  deform  in  the  plane  and  buckle  out- 
of-plane.  Two  sets  of  boundary  conditions  and  rise-to-span  ratios  are 
considered.  In  addition,  results  are  also  presented  for  an  arch  which  first 
buckles  in  its  plane,  sways  to  the  wide,  and  subsequently  buckles  out-of¬ 
plane.  In  Section  5.3.5.,  lateral  buckling  of  a  beam  with  warping  restraint 
is  considered  and  two  examples  are  presented. 

5.2.  Description  of  Problems 

All  the  arches  in  problems  involving  three-dimensional  behavior 
are  assumed  to  have  inextensional  centerlines  and  to  be  fixed  at  the 
boundaries  insofar  as  out-of-plane  motion  is  concerned.  In  certain  of  the 
three-dimensional  problems  selected,  rotations  are  permitted  at  the  supports 
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about  an  axis  perpendicular  to  the  original  plane  of  the  arch  (see  Fig. 
3(a)).  The  two-dimensional  problems  may  involve  either  extensioial  or  in- 
extensional  centerlines  and,  in  addition,  the  arches  may  be  fixed  or  simply- 
supported  at  the  ends.  The  cross  sectional  properties  are  given  in  Table  8. 
All  of  the  arch  members  are  loaded  with  a  concentrated  load  at  the  crown 
(see  Fig.  5(a)). 

In  addition,  some  results  are  presented  for  the  lateral  buciding 
of  an  initially  straight  I-beam  under  a  uniform  dead  load.  Restraint  of 
warping  of  the  cross-sections  is  included  in  the  behavior  of  these  par¬ 
ticular  members.  One  of  the  member  is  a  rolled  steel  section  16  WF  64  and 
the  other  is  a  section  especially  contrived  to  demonstrate  a  particular 
point.  The  cross  section  of  this  special  member  is  shown  in  Fig.  6(b). 

Unless  otherwise  noted,  all  buckling  loads  are  of  the  bifurcation 
type  as  opposed  to  limit  points.  The  following  notation  is  used  in  the 
Figures  and  Tables. 

a  =  non-dimensionalized  load  for  out-of-plane  buckling 
problems,  Q  =  pa2/^^ 

8  =  non- dimensionalized  load  for  in-plane  buckling 
proles,  e.?l2/sl 

JL 

H  =  rise  of  arch 

L  =  span  of  arch 

1^  =  for  a  planar  member,  moment  of  ir  rtia  about  an 
axis  perpendicular  to  plane 

Ij  =  for  a  planar  member,  moment  of  inertia  about  the 
axis  in  the  plane 

J  =  St.-Venart  torsion  constant 
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Cw  =  warping  constant 

e  =  strain  of  centerline  of  member 
c 

5.3.  Numerical  Results 

5.3.1.  Prediction  of  Buckling  Loads 

Data  are  given  in  Table  2  which  indicate  the  rate  of  convergence 
of  the  process  of  predicting  bifurcations.  In  general,  the  change  of  sign 
of  the  determinant  of  the  equations  expressing  the  boundary  conditions  is 
used  to  obtain  an  initial  estimate  of  the  buckling  load.  Then  the  predic¬ 
tion  process  is  implemented  to  "home  in"  on  the  actual  value.  As  may  be 
seen  from  the  successive  values  of  P  and  6P  in  Table  2,  it  is  necessary  to 
apply  the  procedure  taking  advantage  of  orthogonality  between  the  eigen¬ 
vector  and  changes  in  the  prebuckling  configuration  near  the  buckling  load 
in  order  to  guarantee  convergence  (see  Sec.  4,5.).  From  Table  2,  the  case 
of  out-of-plane  buckling  is  seen  to  converge  quite  rapidly  even  though  the 
initial  estimate  of  the  buckling  load  is  in  error  by  a  factor  of  more  than 
three.  This  is  to  be  expected,  since  the  problem  is  essentially  a  clas¬ 
sical  eigenvalue  problem.  That  is,  the  prebuckling  deformations  are  of 
relatively  slight  importance. 

The  last  case  given  in  Table  2  indicates  that  it  is  possible  to 
avoid  the  use  of  the  determinant  involving  the  boundary  conditions  in 
isolating  the  buckling  load.  In  this  particular  case,  an  increment  of 
deflection  was  introduced  and  then  the  prediction  process  implemented  far 
from  the  actual  buckling  load.  Although  the  process  is  seen  to  converge, 
it  is  probably  less  efficient  to  start  the  prediction  process  this  far 
from  the  buckling  load. 
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There  are  some  apparent  minor  discrepancies  in  Table  2.  The 
errors  in  the  coordinate  of  the  load,  as  well  as  the  buckling  load  it¬ 
self,  are  somewhat  dependent  (in  the  fourth  or  fifth  significant  figures) 
on  the  number  of  integration  intervals  as  well  as  the  number  of  cycles  of 
the  Newton-Raphson  technique.  Where  a  direct  comparison  is  made  in 
Table  2  (cases  1  and  2)  the  integration  intervals  and  number  of  cycles  of 
Newton-Raphson  are  the  same. 

5.3.2.  Buckling  Loads  and  Deflections  of  Simply  Supported  Arches 

Results  for  the  buckling  loads  and  deflections  of  a  few  typical 
simply  supported  arches  are  given  in  Table  3.  It  is  seen  that  the  results 
agree  well  with  seom  of  the  previous  analytical  and  experimental  work. 
Figures  7  and  8  show  both  the  prebuckling  and  a  part  of  the  postbuclcling 
curve  for  the  simply  supported  arches.  The  results  plotted  are  for  an 
inextensional  centerline  since  the  effect  of  extension  is  negligible  for 
the  simply  supported  arches  studied  here.  From  Figs.  7  and  8,  it  is  seen 
that  for  H/L  equal  to  0.50,  the  load  carrying  capacity  increases  after 
bifurcation.  This  has  been  observed  experimentally  by  Langhaar,  Boresi. 
and  Carver  (1954)  where,  under  a  concentrated  gravity  load,  the  arch  did 
not  collapse  upon  entering  the  side-sway  buckling  mode.  For  H/L  =  0.25, 
the  load  carrying  capacity  of  the  arch  decreases  rapidly  after  buckling 
(see  Figs.  7(b)  and  8(b)).  This  agrees  with  the  analytical  result  of 
Huddleston  (1968).  Figures  7(b)  and  8(b)  indicate  that  the  method  can  be 
used  to  trace  as  much  of  the  postbuckling  configuration  as  desired. 

The  data  given  in  Table  3  indicate  that  the  stiffness  of  a 
simply  supported  arch,  H/L  =  0.50,  is  slightly  reduced  when  extension  of 
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the  centerline  is  permitted  in  prebuckling  and  postbuckling  behavior. 
However,  the  buckling  load  for  this  arch  is  increased  when  extension  is 
taken  into  account.  This  is  not  a  contradiction  of  Rayleigh’s  theorem 
(1894)  since  bifurcations  from  two  different  prebuckling  configurations 
are  being  compared  and  there  is  no  way  to  assess  the  effect  of  the  internal 
constraint  (ec  =  0).  This  phenomenon  of  a  more  flexible  structure  having 
a  higher  buckling  load  was  reported  by  Masur,  Chang  and  Donnell (1961) . 

In  that  study,  a  gable  frame  with  a  concentrated  load  at  the  peak  was 
analyzed  both  with  and  without  an  inextensible  tie  connecting  the  tops  of 
the  columns.  Removal  of  the  tie  results,  of  course,  in  large  prebuckling 
deformations,  but,  paradoxically,  increases  buckling  load.  Another  in¬ 
stance  of  this  same  phenomenon  occurs  in  another  part  of  the  present  study 
concerning  the  out-of-plane  buckling  of  arches  which  are  either  simply 
supported  or  clamped  in  the  plane.  The  simply  supported  arches  given  signi¬ 
ficantly  higher  buckling  loads  than  the  clamped  ones  for  the  same  H/L  even 
though  they  are  more  flexible  than  the  latter  (see  Fig.  10). 

For  in-plane  buckling  problems,  each  cycle  of  Newton-Raphson  re¬ 
quired  about  one  second  of  computer  time  on  an  IBM  360-75  system.  Usually 
two  additional  cycles  of  Newton-Raphson  sufficed  to  decrease  the  residuals 
to  less  than  0.05  percent  of  their  values  computed  at  the  end  of  the  first 
cycle.  In  the  computations,  only  two  load  increments  were  needed  to  arrive 
at  the  vicinity  of  the  bifurcation  point  for  H/L  =  0.25,  and  three  load 
increments  for  H/L  =  0.50. 
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5.3.3.  Two  Dimensional  Arches  with  Clamped  Ends 

Considerable  analytical  and  experimental  work  has  been  done  on 
shallow  clamped  arches.  Cne  of  the  sample  problems  in  this  study  was  solved 
for  comparison  with  the  experimental  work  of  Gjelsvik  and  Bodner  (1962)  and 
the  analytical  work  of  Schreyer  and  Masur  (1966)  on  shallow  arches  under 
concentrated  loads.  As  may  be  seen  from  Table  3,  the  comparison  with  the 
results  given  by  Schreyer  and  Masur  is  quite  good.  The  agreement  with  the 
experimental  work  of  Gjelsvik  and  Eodner  is  not  as  close,  but  there  are 
uncertainties  in  the  experiments  involving  support  conditions,  modulus  of 
elasticity,  loading  and  dead  weight  of  the  arch.  It  is  appropriate  to 
point  out  that  Gjelsvik  and  Bodner  recorded  the  buckling  load  as  a  maximum 
on  the  experimental  load-deflection  curve  whereas,  the  buckling  load  com¬ 
puted  here  is  of  the  bifurcation  type  and  occurs  after  the  limit  point  (see 
Fig.  9(a))  on  the  load-deflection  curve.  Schreyer  and  Masur  noted  that 
arches  with  certain  rise-to-span  ratios  exhibit  rhis  phenomenon  of  bifurca¬ 
tion  buckling  after  P  falls  off  from  the  value  at  a  limit  point.  As 
expected,  extension  of  the  centerline  is  significant  for  shallow  clamped 
arches,  as  may  be  seen  from  Table  3. 

Results  are  also  presented  for  a  rather  steep  clamped  arch 
(H/L  =  0.25)  which  does  not  buckle  but  rather  maintains  a  symmetrical 
configuration  during  the  loading  process  (see  Fig.  9(b)). 

5.3.4.  Buckling  Loads  and  Displacements  for  Three-Dimensional  Arches 

Table  4  gives  non-dimensionalized  data  for  the  buckling  loads  of 
four  stmple  problems  of  out-of-plane  buckling  of  initially  planar  arches. 
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No  resulto  were  found  in  the  literature  with  which  to  compare  these  results 
directly.  However,  Timoshenko  and  Gere  (1961)  present  some  results  for 
the  out-of-plane  buckling  of  a  uniformly  compressed  arch  which  seem  con¬ 
sistent  with  the  results  obtained  here. 

For  a  given  H/L,  the  simply  supported  arches  have  a  higher  buckling 
load  than  the  clamped  arches,  although  the  clamped  arches  are  initially 
stiffer.  As  may  be  seen  from  Figs.  10  and  11,  all  the  arch  members  examined 
in  this  study  had  reserve  load  carrying  capacity  after  the  buckling  load 
was  reached. 

In  Table  5,  results  are  given  for  an  arch  with  a  section  devised 
so  that  it  first  buckles  in  the  plane  and,  upon  continued  loading,  later 
buckles  out-of-plane.  To  conserve  computer  time,  40  points  on  the  arch  were 
used  in  this  problem  instead  of  100  in  the  numerical  integration  process. 

This  is  the  reason  for  the  slight  discrepancy  between  the  results  presented 
for  this  problem  and  for  the  two-dimensional  problems.  Figure  6(a)  is  a 
schematic  of  what  the  member  cross  section  might  be  in  order  to  have  the 
required  relationships  among  the  three  rigidities. 

5.3.5.  Lateral  Buckling  of  I-Beams 

Results  are  given  in  Table  6  for  the  lateral  buckling  load  of  a 
clamped  I-beam  under  a  uniform  load.  As  may  be  seen  from  Table  6,  the 
result  is  in  excellent  agreement  with  the  previous  work  by  Austin,  Yegiau 
and  Tung  (1957).  The  suppression  technique  is  used  here  to  derive  these 
results.  It  appears  that  the  lateral  buckling  analysis  of  most  rolled  beams 
may  proceed  straightforwardly  as  an  initial  value  problem  without  resorting 
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to  use  of  the  suppression  technique.  Allowable  stresses  and  deflections 
preclude  extremely  long  members  which  give  rise  to  numerical  difficulties. 
When  the  rolled  sections  are  used  as  arches,  however,  the  loads  can  be 
partially  supported  by  normal  forces.  This  makes  possible  a  longer  member 
and  increases  the  effect  of  unwanted  growing  solutions  during  the  numerical 
integration  process. 

Thus,  there  are  cases  in  which  some  technique  like  suppression  is 
required  in  order  to  obtain  accurate  answers,  even  with  double  precision 
arithmetic.  The  numerical  difficulty  arises  when  the  net  effect  of  warping 
restraint  on  the  torsional  stiffness  of  the  whole  member  is  small.  In  this 
case,  the  warping  restraint  is  only  an  edge  effect.  A  long,  slender  member 
is  then  indicated  if  a  computation  is  to  be  carried  out  to  indicate  what 
the  consequences  of  growing  solutions  might  be.  The  section  of  Fig.  6(b) 
was  used  as  a  long  beam  and  the  lateral  buckling  load  sought.  Results  are 
given  in  Table  6  for  the  buckling  load  of  the  member  and  are  given  in 
Table  7  for  a  comparison  of  the  behavior  of  the  solution  versus  the  number 
of  suppressions  used.  As  may  be  seen  from  Table  7,  ten  suppressions  are 
sufficient  to  ensure  satisfaction  of  the  boundary  conditions  while  two  sup¬ 
pressions  lead  to  diverging  approximations. 

Although  results  are  not  given  here,  as  a  matter  of  curiosity, 
the  beginning  of  the  postbuckling  curve  for  lateral  buckling  of  an  I-beam 
under  a  uniform  load  was  computed.  For  the  particular  member,  the  load 
carrying  capacity  dropped,  off  after  buckling.  This  behavior  seems  quite 
reasonable  since  the  late’’  l  buckling  is  accompanied  by  rotation  of  the 
cross  section,  bringing  the  smaller  flexural  igidity  into  prominence. 
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6.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  STUDY 

6.1.  Summary  of  the  Computational  Procedures 

T  a  methods  developed  in  this  study  for  the  analysis  of  buckling 
and  postbuckling  behavior  can  be  summarized  as  follows.  A  new  method  is 
presented  in  Chapters  2,  3,  and  4  for  improving  an  initial  approximation  to 
a  bifurcation  point  on  a  nonlinear  load-def lection  curve.  In  addition,  an 
approximation  to  the  eigenvector  is  generated  and  improved  simultaneously 
with  the  prebuckling  config  iration.  The  initial  stages  of  postbuckling 
are  investigated  by  adding  a  multiple  of  the  eigenvector  to  the  prebuckling 
configuration  at  the  onset  of  buckling.  Subsequent  postbuckling  behavior 
may  be  examined  by  the  application  of  the  standard  Newton-Raphson  procedure 
as  described  in  Chapter  3. 

Tt  >  numerical  methods  introduced  here  for  solving  buckling  and 
postbuckling  problems  involve  two  mod4 Cicat ions  of  the  usual  Newton-Raphson 
technique.  The  first  of  these  modifications  extends  the  Newton-Raphson 
technique  to  the  simultaneous  improvement  of  eigenvalues  and  eigenvectors 
when  there  is  no  difficulty  in  computing  changes  in  the  prebuckling  config¬ 
uration  accurately.  As  indicated  in  Chapter  2,  a  difficulty  occurs,  in 
general,  in  the  vicinity  of  bifurcation  points  where  the  equations  specifying 
changes  in  the  prebuckling  configuration  become  ill-conditioned.  A  second 
modification  of  the  usual  Newton-Raphson  technique  has  been  devised  to  per¬ 
mit  calculation  of  changes  in  the  prebuckling  configuration  in  the  neighbor¬ 
hood  of  a  bifurcation.  In  this  variant  of  the  procedure,  the  orthogonality 
relation  between  the  eigenvector  and  changes  in  the  prebuckling  configura¬ 
tion  plays  an  essential  role. 
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6.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  STUDY 

6.1.  Summary  of  the  Computational  Procedures 

?  a  methods  developed  in  this  study  for  the  analysis  of  buckling 
and  postbuckling  behavior  can  be  summarized  as  follows.  A  new  method  is 
presented  in  Chapters  2,  3,  and  4  for  improving  an  initial  approximation  to 
a  bifurcation  point  on  a  nonlinear  load-deflection  curve.  In  addition,  an 
approximation  to  the  eigenvector  is  generated  and  improved  simultaneously 
with  the  prebuckling  config  iration.  The  initial  stages  of  postbuckling 
are  investigated  by  adding  a  multiple  of  the  eigenvector  to  the  prebuckling 
configuration  at  the  onset  of  buckling.  Subsequent  postbuckling  behavior 
may  be  examined  by  the  application  of  the  standard  Newton-Raphson  procedure 
as  described  in  Chapter  3. 

Tt  •  numerical  methods  introduced  here  for  solving  buckling  and 
postbuckling  problems  involve  two  mod4rications  of  the  usual  Newcou-Raphson 
technique.  The  first  of  these  modifications  extends  the  Newton-Raphson 
technique  to  the  simultaneous  improvement  of  eigenvalues  and  eigenvectors 
when  there  is  no  difficulty  in  computing  changes  in  the  prebuckling  config¬ 
uration  accurately.  As  indicated  in  Chapter  2,  a  difficulty  occurs,  in 
general,  in  the  vicinity  of  bifurcation  points  where  the  equations  specifying 
changes  in  the  prebuckling  configuration  become  ill-conditioned.  A  second 
modification  of  the  usual  Newton-Raphson  technique  has  been  devised  to  per¬ 
mit  calculation  of  changes  in  the  prebuckling  configuration  in  the  neighbor¬ 
hood  of  a  bifurcation.  In  this  variant  of  the  procedure,  the  orthogonality 
relation  between  the  eigenvector  and  changes  in  the  prebuckling  configura¬ 
tion  plays  an  essential  role. 
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The  suppression  technique  or  some  equivalent  scheme  may  be  neces¬ 
sary  when  numerical  integration  procedures  are  used  to  solve  eigenvalue 
problems  of  plate  and  shell  structures.  It  is  well  known  that  *he  differen¬ 
tial  equations  expressing  the  behavior  of  plate  and  shell  structures  have 
edge  effects  as  part  of  their  solution.  A  technique  such  as  the  shooting 
method  would  be  especially  difficult  to  apply  to  such  problems. 

Although  the  numerical  examples  were  chosen  primarily  to  demon¬ 
strate  the  capabilities  of  the  numerical  technique,  seme  interesting  behavior 
of  various  arches  has  been  found.  It  appears  that  in  some  cases  a  more 
flexible  structure  (in  so  far  as  prebuckling  deformations  are  concerned)  may 
have  a  higher  buckling  load.  This  wsj  observed  in  the  in-plane  buckling  of 
an  initially  planar,  simply-supported  arch  under  a  concentrated  load.  When 
extension  of  the  centerline  was  permitted  the  buckling  load  was  higher  than 
its  counterpart  when  extension  was  restrained.  Similarly,  in  the  out-of- 
plane  buckling  of  an  initially  planar  arch,  for  a  given  H/L,  the  simply  sup¬ 
ported  arches  had  a  higher  buckling  load  than  the  clamped  arches.  The  ef¬ 
fectiveness  of  the  numerical  techniques  is  indicated  in  a  particularly 
striking  manner  by  the  somewhat  artificial  problem  of  the  special  arch  member 
discussed  in  Chapter  5  which  buckled  in  its  plane  first  and  subsequently 
out-of-plane.  No  difficulty  was  experienced  in  following  this  complex  load- 
deflection  path. 

6 , 3.  Recommendations  for  Further  Study 

The  proposed  method  may  be  applied  to  many  practical  problems  of 
technical  interest.  Buckling  and  vibrations  of  thin  carved  members  where 
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initial  stresses  are  present  can  be  treated  with  minor  changes  in  the 
computer  codes  develop’d  in  this  study.  In  addition,  nonlinear  stress-strain 
lavs  could  be  admitted  where  the  problem  precludes  significant  unloading. 

The  method  may  also  be  extended  to  eigenvalue  problems  in  plate 
and  shell  type  structures.  The  general  procedure  is  unchanged.  However, 
the  linearized  problems  must  be  treated  by  a  technique  for  approximate 
•solution  of  linear  partial  differential  equations,  rather  than  ordinary 
differential  equations . 

Certain  eigenvalue  problems  in  ryroscopic  motion  may  also  be 
solved,  as  is  obvious  from  Kirchhoff 's  kinetic  analogue  and  the  general 
theory  developed  here  (see  Kirchhoff  (1859)  and  Love  (1927)). 

The  problem  of  deciding  which  boundary  condition  to  relay  when 
generating  the  approximate  eigenvector  needs  -.re  study.  A  sure,  but  some¬ 
what  Inelegant,  solution  to  this  difficulty  is  suggested  in  Section  4.6. 
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TABLE  1.  INITIAL  VALUES  AND  RESIDUALS  FOR  CLArfPED  ARCH 


Quantity 

Homogeneous  Solutions 

Particular 

Solution 

1 

2 

3 

4 

5 

6 

7 

8 

6N1 

1 

0 

0 

0 

0 

0 

0 

0 

W2 

0 

1 

0 

0 

0 

0 

0 

0 

5N3 

0 

0 

1 

0 

0 

0 

0 

0 

6M], 

0 

0 

0 

la 

0 

0 

0 

0 

6MZ 

0 

0 

0 

0 

1 

0 

0 

0 

5M3 

3 

c 

0 

0 

0 

1 

0 

0 

6P 

0 

0 

0 

0 

0 

0 

lb 

0 

Right-Hand- Sides  of 
Linearized  Equations 

0 

0 

0 

0 

0 

0 

0 

R 

Corresponding  initial  incremental  curvatures  are  computed  by  use  of  Eq.  (3.14) 

Jj 

Not  really  a.i  initial  value  since  it  enters  the  computations  at  concentrated 
load  in  middle  of  member  (Eq.  (3.14)) 

I 

I 

| 


i 


TABLE  2.  PREDICTION  OF  BUCKLING  LOADS 
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TABLE  3.  IN-PLANE  BUCKLING  LOADS  OF  ARCHES 


Boundary 

Conditions 

Extension  of 
Centerline 

K/L 

Pa^/EI^ 

d2/L 

Source 

Simply 

Supported 

yes 

.25 

12.981 

.06815 

present 

Simply 

Supported 

no 

.25 

13.006 

.06727 

present 

.'imply 

Supported 

no 

.25 

13.05 

— — 

b 

„ Imply 
Supported 

no 

.25 

13.0 

— — — 

c 

Simply 

Supported 

yes 

.50 

5.8703 

.09762 

present 

Simply 

Supported 

no 

.50 

5.8685 

.09746 

present 

Simply 

Supported 

a 

yes 

.50 

6.54 

— - 

c 

Simply 

Supported 

yes 

.50 

6.15 

- - 

c 

(experiment) 

Simply 

Supported 

a 

yes 

.50 

5.6 

e 

Simply 

Supported 

no 

.50 

5.86 

b 

Clamped 

yes 

.044 

71.866 

.02565 

present 

Clamped 

no 

.044 

77.777 

.02206 

present 

Clamped 

yes 

.044 

72.2 

— 

f 

Clamped 

yes 

.044 

63.7 

g 

(experiment ) 

Extension  of  the  centerline  was  permitted  in  deriving  the  prebuckling 
configuration,  but  not  in  the  eigenvector. 

^Schmidt  (1969) 


CHuddleston  (1967) 

^Langhaar,  Boresi,  and  Carver  (1954) 
eChen  and  Boresi  (1961) 

^Schreyer  and  Masur  (1966) 

8G.ielsvik  and  Bcdner  (1962) 


Clamped 


Simply-Supported 


Pa2/EI2  GJ  d2/a  Pa2/EI2  GJ  d2/a 


3.453  0.0007857  3.952  0.001389 

0.6684  0.0003262  0.7701  0.0C06080 


TABLE  5.  BUCKLING  LOADS  AND  DEFLECTIONS  FOR  A  SIMPLY- 
SUPPORTED  ARCH  WHICH  FIRST  BUCKLES  IN-PLANE 
AND  UPON  INCREASED  LOADING  SUCKLES  OUT-OF¬ 
PLANE,  H/L  »  0.25,  e  =  0 

c 


In-plane  Bucklin 


Pa2/EI, 


Subsequent  Out-of-Plane  Bucklin 


d2/L  Pa  /EI^ 


13.04c 


0.06648  12,70 


0.07851 


0.3562 


aThis  differs  from  the  results  for  the  in-plane  buckling  of  other 
two  dimensional  arches  because  fewer  points  were  used  here  in  the 
numerical  integration  process. 


TABLE  6.  LATERAL  BUCKLING  LOADS  OF  UNIFORMLY  LOADED,  CLAMPED  I-BEAMS 
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TABLE  8.  MEMBER  SECTION  PROPERTIES 


Out-of-plane  buckling  I. 

(12WF31) 

In-plane  buckling  of  1^ 

simply  supported  arches 

In-plane  buckling  of  1^ 

clamped  arches 


238.4  in4,  I2  =  19.8  in4,  J  =  .5065  in4 

18.0  in4,  Area  =6.0  in2 

.5493  x  10"3  in4,  Area  =  .1875  in2 


FIG.  1.  GLOBAL  AND  LOCAL  COORDINATE  SYSTEMS 


A,  Deflection 


FIG.  2.  QUALITATIVE  FORCE-DEFLECTION  CURVE 


p 


(b)  Schematic  of  Anti-Symmetrical 
In-plane  Buckling  Mode 


FIG.  5.  TYPICAL  IN-PLANE  BEHAVIOR  OF 
SIMPLY  SUPPORTED  ARCH 


(a)  Schematic  of  Cross-Section  for 
Member  Which  Buckles  In-Plane 
and  Then  Out-of-Plane 


(b)  Schematic  cf  Cross-Section  for  Special 
Member  in  Lateral  Buckling  Study 


FIG.  6.  SPECIAL  CROSS  SECTIONS  OF  MEMBERS  USED  IN  THE  ANALYSIS 


7.  LOAD  VERSUS  VERTICAL  DEFLECTION  AT  CROWN,  IN-PLANE 
BUCKLING  OF  SIMPLY  SUPPORTED  ARCHES,  e  =  0 


0 


.02 


.04 


.06 


.08 


(b)  H/L  -  .25  d3/L 


FIG.  8.  LOAD  VF  'US  HORIZONTAL  DEFLECTION  AT  CROWN,  IN-PLATE 
BUCKLING  OF  SIMPLY  SUPPORTED  ARCHES,  z^  *  0 


I 


T 


T 


(a)  H/L  =  .25 


(b)  H/L  =  .5 


10.  LOAD  VERSUS  VERTICAL  DEFLECTION  AT  CROWN,  OUT-OF-PLANE 
BUCKLING  OF  CLAMPED  AND  SIMPLY  SUPPORTED  ARCHES,  s  =  0 


Suoported 


FIG.  11.  LOAD  VERSUS  OUT-OF-PLANE  DEFLECTION  AT  CROWN,  OUT-OF-FLAME 
BUCKLING  OF  CLAMPED  AND  SIMPLY  SUPPORTED  ARCHES.  *  =  0 
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APPENDIX  A 

SOLVABILITY  OF  THE  BASIC  EQUATIONS  OF  THE  METHOD 
A.l.  Case  of  a  Single  Root 

Consider  the  method  proposed  in  Chapter  2  as  applied  to  the 
determination  of  the  bifurcation  point  and  corresponding  eigenvector  ce  the 
algebraic  system 


A  X  =  X  B  X  (A.l) 

where  for  purposes  of  this  discussion,  A,  B,  and  X  correspond  to  the  onset 
of  buckling.  At  buckling,  both  A  and  B  are  assumed  to  be  self-adjoint  and 
B  is  taken  to  be  positive  definite.  The  side  condition,  corresponding  to 
Eq.  (2.6)  is  taken  as 

XT  B  6X  =  0  (A. 2) 


At  the  buckling  point,  the  coefficient  matrix  given  by  tin’  left-hand-side 
of  Eq,  (2.3)  and  Eq.  (A  2)  is 


C  = 


A  -  XcrB 


T 

f-  x£  B 


L_ 


I 


B  x. 
J. 

0 


(A  >3) 


where  X  is  the  buckling  load  and  is  the  corresponding  eigenvector. 

The  basic  method  will  fail  if  the  coefficient  matrix  C  of  Eq,  (A, 3), 
used  in  the  compute tion  of  the  increments  of  an  approximate  eigenvalue  and 
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eigenvector,  is  singular.  It  is  expected  that  if  this  occurs,  the  singu¬ 
larity  will  exist  at  exactly  the  prebuckling  configuration  given  by  A,  B, 
and  Acr*  If  the  order  of  the  original  problem  is  of  order  n,  then  C  in 
Eq.  A.  3  is  a  symmetrical  matrix  of  order  n+1. 

The  matrix  C  in  Eq.  A. 3  will  now  be  shown  to  be  nonsingular  by 
a  consideration  of  the  eigenvalues  of  the  auxiliary  system 


C  y  =  A  D  y 


(A. 4) 


where 


(A. 5) 


It  may  be  verified  by  direct  substitution  that  the  eigenvectors  y  , 
(m  «  1,  .  n+1)  of  the  system  given  by  Eq.  A. 4  are 


(k  =  2,  ...... n)  where  the  x^  and  x^  are 


eigenvectors  of  Eq.  (A.l).  The  corresponding  eigenvalues  A  of  Eq.  (A. 4) 
are  1,  +1,  and  (Ak  -  Afr).  The  eigenvectors  of  Eq.  (A.l)  are  found  by  con¬ 
sidering  A  and  B  constant  at  the  prebuckling  configuration  corresponding  to 
the  onset  of  buckling,  and  are  assumed  to  be  normalized  with  respect  to  B. 

It  is  not  difficult  to  show  that  the  determinant  of  C  is  equal  to 
the  product  of  the  A's  multiplied  by  det  (D) .  Since  the  latter  is  equal  to 
det  (B)  which  is  positive,  then  det  (C)  is  nonzero  provided  none  of  the  A 
are  zero.  Only  in  the  case  of  a  multiple  root  can  a  A  be  zero.  Thus,  if 
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there  are  no  multiple  eigenvalues  of  the  original  system  given  by  Eq.  (A.l), 
the  basic  method  proposed  encounters  no  numerical  difficulties  associated 
with  a  singularity  of  C. 

A. 2.  Case  of  a  Double  Root 

The  existence  of  a  double  root  of  Eq.  (A.l)  (say  X  =  A,r)  implies 

cr  K 

that  the  matrix  C  in  Eq.  (A. 3)  is  singular  at  the  bifurcation  point,  xhis 
singularity  may  be  removed  by  the  following  computational  sequence.  Two 
independent  eigenvectors  are  generated  by  specifying  two  side  conditions  £or 
each  eigenvector.  The  two  eigenvectors  are  denoted  here  by  and  x^  and 
their  increments  by  6x^  and  6x^.  The  side  conditions  for  6x^  are 

x*  B  6x^  =  0  ,  xi  B  Sx^  “  0  (A. 6) 

and  the  side  conditions  for  dx^  are 

x*  B  dx,,  =  0  ,  x^  B  dx  =  0  (A. 7) 

K.  K  JL  is. 

The  specification  of  the  two  side  conditions  results  in  the  fol¬ 
lowing  coefficient  matrix  for  the  equations  determining  the  incremental 
changes  in  the  two  eigenvectors 


80 


where  D  is  given  by  Eq.  (A. 5),  C  is  given  by  Eq.  (A. 3),  and  yR  =  \  -g- 

The  coefficient  matrix  C  has  one  more  row  than  column,  bat  as  indicated  by 

Koiter  (1945) ,  the  equations  which  give  rise  to  C  are  all  valid  at  the 

bifurcation  point.  An  independent  set  of  equations  with  a  nonsingular 

—  — T 

coefficient  matrix  may  be  derived  by  premultiplying  C  by  C  .  The  result 
of  this  multiplication,  which  amounts  to  an  application  of  a  least  squares 
technique,  is 


?C  =  0T  0  +  D  ,KyJ  D  (A.9) 

The  object  is  to  show  that  the  coefficient  matrix  in  Eq.  (A.9)  is 

nonsingular.  The  eigenvector  yK  corresponds  to  a  zero  eigenvalue  of  the 

matrix  C  of  Eq.  (A. 4).  As  shown  in  Section  A.l,  the  remaining  eigenvalues 

T 

of  C  are  nonzero.  The  matrix  C  C  in  Eq.  (A.9)  has  the  same  eigenvectors  as 

T 

C.  It  follows  that  the  eigenvalues  of  C  C  are  the  squares  of  those  of  C. 

T 

Now  consider  the  matrix  G  =  in  Eq.  (A.9).  Direct  substitution  yields 

the  result 


G  yK  =  1  D  yR  (A. 10) 

From  Eq.  (A. 10)  it  may  be  seen  that  the  eigenvector  is  also  an  eigen¬ 
vector  of  G  and  the  corresponding  eigenvalue  is  unity.  The  matrix  G  is 
constructed  in  such  a  way  that  its  remaining  eigenvalues  are  zero  since  it 
is  a  symmetric  matrix  of  rank  one.  The  remaining  eigenvectors  of  G  may 
therefore  be  taken  the  same  as  those  of  C. 
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T  T 

Thus  both  matrices  C  C  and  Dy^y^D  in  Eq.  (A. 9)  have  the  same 

eigenvectors.  The  eigenvalue  of  the  sum  of  two  matrices  having  the  same 

eigenvectors  is  merely  the  sum  of  the  eigenvalues  of  the  individual  matrices. 

— T—  T 

It  follows  that  the  eigenvalues  of  C  C  are  those  of  C  C  except  for  the  zero 

T 

eigenvalue  which  becomes  +  1  (from  the  matrix  Dy^y^D) .  Since  all  the 
— T— 

eigenvalues  of  C  C  are  nonzero,  it  is  nonsingular  and  the  method  proceeds 


without  difficulty. 


APPENDIX  B 


SOLVABILITY  OF  THE  EQUATIONS  USED  IN  DETERMINING  ACCURATE  CHANGES 
IN  THE  PREBUCKLING  CONFIGURATION  NEAR  A  BIFURCATION  POINT 

The  linearized  operator  used  to  compute  changes  in  the  prebuckling 
configuration  becomes  singular  at  bifurcation  points,  as  has  been  noted  by 
Thurston  (1968).  This  singular  operator  is  denoted  here  by  D  where 

D  =  A  -  X  B  (B.l) 

cr 

The  discussion  here  will  be  limited  to  the  algebraic  eigenvalue  problem  so 
that  A  and  B  are  matrices  which  define  the  prebuckling  configuration  at 
the  onset  of  ’-uckling  and,  X£r  is  the  buckling  load.  The  matrices  A  and  B 
are  assumed  to  be  self-adjoint  and  B  is  taken  to  be  positive  definite. 

A  technique  has  been  discussed  in  Chapter  2  for  removing  the 
singularity  from  D.  It  is  the  object  of  this  Appendix  to  show  that  the 
resulting  coefficient  matrix  is  indeed  nonsingular.  As  indicated  in 
Chapters  2,  3,  and  4,  a  side  condition  is  appended  to  the  basic  system. 

This  side  condition  specifies  that  changes  in  the  prebuckling  configuration 
are  orthogonal  to  the  eigenvector  and  may  be  expressed  formally  as 

xj  B  y  =  0  (B.2) 

where  x^  is  the  eigenvector  corresponding  to  the  singularity  of  B  and  y  is 
the  change  in  the  prebuckling  configuration.  This  side  condition  leads 

'V, 

to  a  new  coefficient  matrix  D  given  by 
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2f  = 


1 

A  -  X  B 

_ £E__  ! 

T  _  I 
X1B 


(B.3) 


which  has  one  more  row  than  column.  As  mentioned  in  Chapter  2,  all  these 

'V 

equations  giving  rise  to  D  are  valid  at  Che  bifurcation  point. 

A  consistent  set  of  equations  with  a  nonsingular  coefficient  matrix 
is  derived  using  the  least  squares  technique: 


DTD  =  (A  -  X  B)T  (A  -  X  B)  +  Bx,x:F  B 
cr  cr  il 


(B.4) 


The  matrix  given  in  Eq.  (B.4)  may  be  shown  to  be  nonsingular  by  an  argument 
exactly  parallel  to  that  given  in  Appendix  A,  Section  A. 2  for  the  case  of  a 
double  eigenvalue. 
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APPENDIX  C 

ENSURING  ORTHONORMALITY  OF  THE  DIRECTION  COSINES 

The  particular  technique  used  in  this  study  for  handling  the 
geometry  treats  each  of  the  nine  direction  cosines  as  an  independent  quan¬ 
tity  during  certain  stages  of  the  numerical  computations.  Since  the 
direction  cosines  are  required  to  form  an  orthonormal  set,  it  is  necessary 
to  enforce  this  constraint  in  some  manner.  The  method  for  ensuring  ortho- 
r.ormality  of  the  direction  cosines  is  outlined  below. 

Orthonorm&lity  of  a  set  of  direction  cosines  U  requires  that 

U  UT  =  I  (D.l) 

where  I  is  the  identity  matrix.  Substitution  of  an  approximately  ortho¬ 
normal  set  of  direction  cosines,  U  ,  into  Eq.  (D.l)  yields 

3 

!J  UT  =  I  +  e  S  (D.2) 

a  a 

where  S  is  a  symmetric  error  matrix  whose  individual  elements  are  presumed 
to  be  of  order  urity  and  e  is  small.  A  correction  matrix  C  is  introduced 
such  that 


U  +  C  =  U 

a 


(D.3) 


The  matrix  C  is,  of  course,  not  unique.  A  convenient  choice  is 


C  =  1/2  e  S  U 


(D.4) 
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By  direct  substitution,  it  may  be  shown  that  Eqs.  (D.3)  and  (D.A)  satisfy 

2 

Eq.  (D.2)  to  terms  of  order  e  .  Since  the  quantity  U  in  Eq.  (D.4)  is  not 

known,  U  is  used  as  a  first  approximation  to  U.  Equation  (D.4)  becomes 
Si 

C  *  1/2  e  S  U  (D.5) 

a 

Equation  (D.3)  may  be  used  to  describe  an  iterative  process  where 
U  is  interpreted  as  the  latest  approximation  and  U  as  the  previous  approxi- 
mation  to  the  required  orthonormal  set.  Substitution  of  Eq.  (D.5)  into 
Eq.  (D.3)  and  rearrangement  yields 

U  =  (I  -  1/2  e  S)  U  (D.6) 

a 

At  a  particular  iterative  step,  the  value  of  U  computed  in  Eq.  (D.6)  becomes 
U  for  the  next  step.  When  the  coefficient  e  becomes  small  enough,  the 
correction  process  is  terminated.  This  correction  process  is  necessary  at 
each  integration  point  along  the  member. 

The  correction  process  discussed  above  results  in  a  new  set  of 
direction  cosines  which  is  n  t  derivable  from  the  first  derivatives,  i.e., 

\  d i.  . 

/  C-p1)  ds  4  l±y  (i,  j  =  1,  2,  3)  (D.7) 

The  following  computational  scheme  was  devised  to  ensure  that  Eq.  (D.7) 
is  satisfied.  The  corrected  direction  cosines  are  substituted  into  Eqs. 
(3.3)  and  new  first  derivatives  computed.  A  quadrature  of  these  first 
derivatives  yields  new  direction  cosines  consistent  with  Eq.  (D.7). 
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Residuals  are  computed  from  Eqs.  f3,3)  using  the  direction  cosines  from 

the  quadrature.  The  residuals  are  then  used  to  drive  the  linearized 

geometric  equations  of  Chapter  3. 

This  technique  has  been  implemented  as  part  of  the  solution  of 

the  geometric  equations  of  Chapter  3.  Before  this  technique  was  devised, 

it  was  not  possible  to  achieve  global  equilibrium  even  though  the  residuals 

in  the  differential  equations  were  small. 

The  effect  of  the  technique  is  to  transfer  the  residuals  in 

Eq.  (D.2)  back  to  the  geometric  differential  equations,  those  of  Chapter  3- 

That  is,  a  residual  of  order  e  in  the  algebraic  equations  results  in  a 

residual  cf  order  e  in  the  differential  equations.  The  application  of 

Newton-Raphson  to  the  differential  equations  gives  rise  to  changes  of 

2 

order  e  in  the  direction  nosines  leading  to  new  residuals  of  order  e  in 


the  direction  cosines. 
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